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Introduction. While considerable progress has been made in recent years in various 
stability problems of hydrodynamics (see C. C. Lin [1] and the bibliography contained 
in the reference) and in the discussion of fully developed turbulence (see G. K. Batchelor 
[2]), the problems associated with transition from laminar to turbulent motion are in a 
less satisfactory state. To be sure the phenomenological theory of H. Emmons [3, 4] 
shows considerable promise, but it replaces explicit dependence on the hydrodynamical 
equations by general stochastic and probabilistic considerations. 

Unfortunately, the usual laminar stability theory of the Orr-Sommerfeld equation 
seems incapable of reasonable extension to the nonlinear regime. In addition, even the 
linear stability analysis is very difficult and requires delicate mathematical consideration. 

In contrast to this, the techniques used by E. Hopf [5] to establish the existence of 
a weak solution for all time of the Navier-Stokes equations in a bounded region are in 
principle simple and constructive, and it therefore seems reasonable to investigate their 
practicality for, first, the investigation of the stability problems of hydrodynamics and, 
second, any insight they might be capable of providing for the transistion problem. 

In essence, Hopf’s method consists of proving the existence of a complete ortho- 
normal set of functions of the space variables which have divergence zero and which 
vanish outside sets of compact support. A solution to the Navier-Stokes equations is 
then sought as an orthogonal series expansion with unknown coefficients of time. This, 
when substituted into the Navier-Stokes equations and sampled by all members of the 
orthogonal set, leads to an infinite set of ordinary nonlinear differential equations in 
the infinite number of unknown time dependent coefficients. The existence of a solution 
of this system is then established by a process of successive approximations involving 
a square truncation of the system. The idea of the weak solution and Hopf’s treatment 
of the boundary conditions closely parallel those of the direct methods of the calculus 
of variation. 

Certainly, if for a given problem, a sufficiently judicious choice of a complete ortho- 
normal set of functions in the space variables could be made so that the main remaining 
features of the Navier-Stokes equations were to be contained in a relatively few non- 
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linear ordinary differential equations, the possibility would then exist that some insight 
would be gained into the mechanism of transition. In fact, it might then be possible to 
verify the appealing intuitive picture of this phenomenon as conjectured by Hopf [6] 
in his model. In brief, he pictures the transition from laminar to turbulent motion as a 
bifurcation phenomenon in phase space in which the number of dimensions of the mani- 
fold of central motions changes from zero (the laminar regime where all perturbations 
die out) by integral values to a large number of dimensions (possibly an infinite number) 
corresponding to fully developed turbulence. 

On the other hand, it could happen that all reasonable and tractable choices for this 
family could, for a given problem, require the consideration of an enormous number of 
equations to represent the Navier-Stokes equations adequately, and it might not be 
practical to replace an existing asymptotic treatment based on the Orr-Sommerfeld 
equations by one based on the ideas of Hopf. That is, the situation could be intrinsically 
like Poincaré’s famous example of a convergent series requiring a million terms for ten 
percent accuracy and only two or three for the asymptotic series. 

In this paper we will demonstrate that the qualitative and quantitative features of 
the instability of plane Poiseuille flow can be obtained from low order truncations of a 
suitably chosen Hopf type of expansion. The qualitative features are already present 
in the approximate curve of neutral stability obtained from an eight-by-eight truncation 
of the infinite system of equations satisfied by the coefficients of the expansion even 
though the minimum critical Reynolds number given by this crude approximate, while 
of the right order of magnitude, is approximately only half that computed by Lin. In 
contrast, the incomplete results presently available for the twenty-by-twenty truncation 
show excellent agreement with the results of Lin with respect to the minimum critical 
Reynolds number. In any event, our methods are elementary and involve only trigono- 
metric and hyperbolic functions and are limited in accuracy only by machine limitations 
as to the size of matrices which can be currently treated. Moreover, as will be shown in 
Sec. II, all solutions of the Orr-Sommerfeld equations are expressible in terms of the 
solution derived here so that the difference between our approximate neutral curve and 
the correct one can be attributed to truncation error. 

We will limit consideration here to the presentation of the technique for those per- 
turbations whose stream functions are even in y, the running variable perpendicular 
to the plate. We are unaware of any detailed calculations for those perturbations whose 
stream functions are odd in y, but certainly they are equally accessible by our methods 
should any doubt really arise as to whether they are all actually stable, as is generally 
supposed. 

The mathematical justification of our techniques, as well as numerical results of 
higher order approximations, will be presented elsewhere. In the main the proofs consist 
of extensions and refinements of several known techniques, some of which have been 
treated for the first time by one of us [7]. Hopf’s original proof is not adequate, since 
Poiseuille flow is not confined to a bounded region. 

In view of the long and often controversial history of this particular problem, we 
feel fortunate to be able to present additional evidence for the instability of plane 
Poiseuille flow by elementary means. Moreover, unlike methods based on energy balance, 
even our eight-by-eight approximations give results of the right order of magnitude. 
Thus we believe that we have made considerable progress in furnishing a simple alternate 
approach to problems of linear instability, the need for which has been stressed recently 
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by Lin [1]. Finally, the qualitative success that we have achieved offers considerable 
promise that consideration of the nonlinear terms would provide insight into the processes 
involved in transition to turbulence. This particular problem may, however, be unsuited 
for the study of realistic nonlinear behavior since, like other workers in the field, we 
have assumed periodicity along the plates. 

As a final introductory remark, we would like to express our sincere thanks to Dr. I. 
Tarnove who carried out the calculations which will be presented here. We would also 
like to extend our thanks to our colleagues, Professors G. Carrier, F. Clauser, W. Hay, 
and L. Lees, and Drs. J. Sellars and G. Solomon, who have encouraged us in this work 
and who have made many constructive suggestions. In particular, the first author is 
indebted to E. Hopf for his hospitality during that week-end in April 1955 when this 
approach was still incompletely formulated (see reference [8]). 

Section I. Formal transformation of the problem into an infinite system of ordinary 
differential equations. The Navier-Stokes equations for the two-dimensional motion of 
an incompressible viscous fluid are written in the familiar form, 


u, + uu, + vu, + p, = vdu, 
v, + uw, + vv, + p, = vdv, (1) 
u, tv, = 0. 


Here u(x, y, t) and v(2, y, t) are, respectively, the z- and y-components of the velocity 
of the fluid particle which, at time ¢, is at the point (z, y). The pressure at this point 
and instant is denoted by p(z, y, t), while v represents the coefficient of viscosity divided 
by the density and A represents the Laplace operator. 

These Navier-Stokes equations are satisfied by the following values for u, v, and p, 
representing so-called plane Poiseuille motion. 


u=1-y’, v= 0, p = —2vx + const. (2) 
We are primarily interested in this solution for -— ~ < 2 < + o and for 
— 1 < y < + 1. It is seen that this solution satisfies the condition requiring that the 
particles of fluid at the boundary should be at rest, namely 
u=v=0 for y= +1. (3) 
The problems connected with the stability of this Poiseuille motion require that we 
consider nearby motions also satisfying the boundary conditions (3). We accordingly 
introduce the transformation, 
u=1—y'+w’, v=v’, p = —2vx + const + p’, (4) 
where for the problem at hand, wu’, v’ and p’ will be thought of as small. In these primed 
variables, the Eqs. (1) appear in the form, 


ul + (1 — yuh — Qyro’ + pt — vAu’ + w’ul + 0'u) = 0. (5a) 
ve + (1 — ol + pl — vdv’ + wu’) + o'r) = 0. (5b) 
u; +v, = 0. (5c) 


It will be useful to eliminate p’ from (5a) and (5b) by differentiating the first with 
respect to y, the second with respect to x, and then forming the difference and simplifying 
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with the aid of (5c). We thus find that 
(uy —v’), +0 — Wu, — of), — 2’ — vA(u) — v!) : 
(6) 
+ u'(ui — vv), + v’'(uh — v}), = 0. 
Equation (5c) permits the introduction of a “stream function,” y¥(x, y, ¢) such that 


u’ = y, andv’ = — y,, and hence (6) may be written in the form, 


a(Ay, Wy) =. (7) 


Ay, + (1 — y)Ay. + 2¥, — vA’y + = 
Or, Y) 


The linear stability theory is based on the equation obtained by neglecting the 
second order terms in (7), namely those which are compressed in the expression 0(Ay, 
y)/a(x, y). We also start out by considering only stability relative to sinusoidal dis- 
turbances of a specified period 27/a in x. Hence, we set y = e***¢(y, ¢) in the linearized 


equations (7) and then divide by e***. The result is 


(6,, —a¢), + (1 —- Y (by — ad)ia + Ziad = v(byywy — 2a°$,y + ad). (8) 


We attempt the solution of (8) with the appropriate boundary conditions derived 
from (3) and the subsequent definitions of y and ¢, namely 


(o, t) = ¢,(0, t) = 0. (9) 
¢ u 


We attempt this solution in the form of an infinite series of the form 


oly, t) = b a,,(t)b,(y) : (10) 
n=1 
where ¢,(y), $2(y), +++ are associated with certain eigenvalues X, , 2 , «++, respectively, 
satisfy the initial conditions (9), and also satisfy the differential equations 
dé > d’o ; 2 dbp 
- = 2a — —- adn = r.( 2%, a a ee (11) 
ay dy dy 


It follows from (11 and (9) that, if Aw * Nas then Din and ¢, are orthogonal to each 
other in the sense that the integral on the left of the next formula (12) is zero. Even if 
X,, could be equal to A, (with m not equal to n), which will turn out to be impossible, 


the Schmidt process of orthogonalization of linearly independent functions permits us 


A 


integral always vanishes when m + n. Finally, assuming hereby 


to assume that the 
that none of the ¢,’s is identically zero, we may introduce normalizing factors, so that, 


without exception, we may write 


[ ladn(yo,.(y) 4 o,,(ydi(y)] dy = Sbnn « (12) 

For purposes oi abbreviation we introduce the notation 
[haw + 1@s'Wlay = 1,9) (13) 
so that (12), for example, can be written in the form (¢,, , ¢.) = Sma - Assuming that 


such an expansion as (10) is possible, that this series can be differentiated term by term, 


and that the differentiated series is uniformly convergent, a simple calculation based 
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on (12) and (13) shows that 


(¢, dn) — a,,(t) = an ad (14) 
We also have a Bessel inequality 
Da < @, 9%), (15) 
k=l 
which establishes the convergence of >-?., a? . Conversely, if a, , a2 , --* were given 


arbitrarily subject to the condition that this series converges, we can establish by known 

and straight-forward methods [9] the existence of a function ¢ vanishing at y = + 1 and 

satisfying (14). In order that ¢ should also satisfy (8), we seek, in a purely formal manner, 

a system of differential equations of infinite order to be satisfied by a,(t), aa(t), +--+ . 
Formally differentiating (10) and substituting in (8), it is found that 


> al(t)[o’(y) — a°o,(y)] + ia(l — y*) >> a,(t)[bh’(y) — 0°¢,(y)] 
n=1 n=l 


(16) 


os) 


+ ia D2 an(idna(y) =v Do an()Anla’d.(y) — o2’(y)]- 
n=1 n=1 

In obtaining this result, the terms on the right-hand side are reduced with the help of 

(11). Integration by parts shows that 


atl 

[ balw)la*ba(y) — 42'Q)] dy = nm 4x) = Snes 

J-1 

inasmuch as ¢,,(-+ 1) = 0 in accordance with the boundary conditions. Hence, on 
multiplying (16) by ¢,,(y) and integrating from — 1 to + 1, we find that 


a+ 


—an() + ia Dat) | (= Pbnlloe’y) — a%a(y)] dy 


J-1 


atl 
+ 2ia >> a,(t) o,(y)b,(y) dy = vr,,a,,(0), m = 1,2,3,- 
1 


n=1 d= 


This equation may be written in the form, 


de = 
a + Wr,.Omn = ia] — i ¥ ena, | , m=1,2,3,-°-, (17) 
( n=l 
where 
a+1 ” 
Cnn = | U'dnly)[adn(y) — on’(y)] dy + 2 [ m(Y)oa(y) dy. (18) 
v—1 = 


It is now necessary to examine a bit more explicitly the eigenvalues \, , A. , As, *** 5 
and the associated eigenfunctions, ¢,(y), ¢2(y), ¢;(y), --* , Which satisfy (11) and vanish 
together with their first derivatives at y = + 1. Equation (11) has only constant co- 
efficients. Hence its general solution is readily found by elementary methods and is 
found to be 


¢(y) = A, ecoshay + A, sinh ay + A; cos(A — a’)'?y + Aysin (A —a’)'”’y, (19) 


where A, , A,, A; , and A, are constants of integration. In the sequel we shall system- 
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atically use the abbreviation 
u=(\—a’)'” Un = (A, — a’)'”?. (20) 


, 


To determine the A’s so that the boundary conditions ¢(+ 1) = ¢’(+ 1) = Oare satisfied, 
we are led to the equations, 


| 


A, cosha + A,sinha + A; cosu + A,sinu = 0, 
A, cosha — A,zsinha + A; cosu — A,sinu = 0, 


0, 


I 


A, asinha + A,a cosha — A; usinu + Ayu cosu 


—A,asinha + A,a cosha + A,;usinu + A,ucosu = 0. 


By forming half the sum and difference of the first pair of these equations and then 
doing likewise for the second pair, we see that this system is equivalent to the following: 


A, cosh a + A; cosu = (, 
A, asinh a — A;,usinu = 0, 
A, sina + A,sinu =0, 
A, a cosh a + A,ucosu = 0. 


The.determinant of this latter system of homogeneous equations in the A’s is easily 
seen to have the value 


fu cosh a sin u + a sinha cos uj[u sinh a cos u — @ cosha sin u]). 


Hence the system is consistent (leading to solutions of (11) that are not identically 
zero and at the same time satisfy the boundary conditions), if and only if 


u cosh asin u + asinh cosu = 0 (21) 

or 
usinh a cosu — a coshasinu = 0. (22) 
It is readily seen that these equations, (21) and (22), cannot be satisfied simul- 
taneously. If (21) is satisfied, A, = A, = 0 and the ¢ given by (19) is an even function. 
If (22) is satisfied, A, = A, = 0 and ¢ is an odd function. This suggests the possibility 
of getting more rapid results by considering separately the even disturbances ¢(y) = 
¢@(— y) and the odd disturbances ¢(y) = — ¢(— y). Any function can be expressed as 


the sum of an odd function and an even function, and hence there is no loss of generality 
in this approach, provided that there is no coupling between the odd and even disturb- 
ances. Since the coefficients of the linear equation (8) are all even, and since only even- 
order derivatives occur with respect to y, it is clear that such a coupling is indeed not 
present in the linearized theory. This fact is also confirmed by formula (18) for the 
“coupling coefficient”’ c,,, , Which is clearly zero if ¢,,(y) and ¢,(y) have opposite parity. 

In view of these considerations we limit ourselves hereafter in this paper to even 
disturbances. Accordingly, we ignore solutions of (22) and denote by u, , uz, Us, --* , the 
successive solutions of (21) and by \, , A», As, «+ the corresponding eigenvalues given 
by (20). We continue to regard a, of course, as fixed. The following may be said about 
the successive solutions of (21). 
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In the first place, Eq. (21) can be written in the form 
utanu = 8, where 8 = —atanha (23) 


and we proceed to get an explicit approximate expression for the solutions of (23) valid 
at least for small values of 8, and hence of a. Evidently (23) is satisfied when 8 = 0 
and u = kr, k = 0, + 1, + 2, --- . By the implicit function theorem, we can solve for 
u as a function of 8, when k = + 1, + 2,---. Wehave 

du u d*u ue + Bu 


.1 B ee 6 en . we _ uw 4 
u— kr = tan s . a of +B+ 3 . dB = 2(8 + 1) Cu +B+B 





-* 


2 


d’u , ; 
“ —2(8 + lu 
dg 


{itht 48+ 7 = SAB! 48+ Wel £83) = WA) 
L (a? + B+ B) 


Hence the first four terms, in the Taylor series development of u as a function of 8, may 





be written down at once, 


2 3 
8 B B 
u=ke+>—-—-paot2mat::: 
kr k T k TT 
and we may estimate the error in retaining only one, two, or three terms with the help 
of the remainder theorem for the Taylor series, together with the above formulas for the 
first three derivatives. In terms of a, then, the kth u is given approximately as follows: 


aod atanha a’ tanh’ a@ a’ tanh’ a 
’ “ kr k’x® k’r® F 
which shows that u, is asymptotically undistinguishable from kz. 
With u, (and hence \, = uz + a’) thus characterized, we define a corresponding 
eigenfunction, 
o*(y) = cosuy — (cos u, cosh ay)/(cosh a), (24) 
by taking A, = 1 arbitrarily in (19) and then determining A, = — (cos u,)/(cosh a) 


so that ¢¢(1) = 0. Of course, A, = A, = 0 as already stated, in the even case to which 
we are restricting attention. The other boundary conditions, viz., ¢¥’(+1) = 0, are seen 
to be satisfied in virtue of the fact that u, , bv definition, satisfies (21). These ¢*’s are 
not quite the same as the ¢’s previously introduced because they are not normalized. 
In order to carry out the normalization, we need to calculate the integral expression 
for (¢* , ¢*) = u,, . This may be done in several obvious ways, using (21), perhaps to 


simplify the result. We thus find that 
Um = A, [1 + (sin 2u,,)/2u,,]. (25) 
Then we also have that 
(% ,O%) = Srnbtn « (26) 
Since ¢, is proportional to ¢* and since (¢,, , d.) = 4Snn , it is evident that 


on = Um On - (27) 
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On account of the awkwardness of the expression for the normalizing factor, it is some- 
times desirable to expand directly in terms of the unnormalized ¢*’s. Thus, instead of 
(10) we may sometimes consider an equivalent expression, 


=) 


o(y, t) = >) ak(D¢*(y), (28) 
n=! 
where 
ax = un (d, O%) = Un “(G, om) = bn Om - (29) 


Correspondingly, we sometimes write (17) in the form 
day. * : * = ‘ 
“dt + vrv,,0% = 1a — Om + ) & Cats ? (30) 
n=1 


where 


2 


re a Se (31) 


It follows from (27) and (31) that (18) can be written in the form 
atl 
ch. = un | y'ok(y)la*or(y) — o8’"(y)] dy 
v—] (32) 
al 
+ 2u,,’ | Pr(yo>r(y) dy. 
-1 
Finally, the integrals appearing on the right-hand side of (32) can be evaluated by 
elementary methods using (21) and (24). We exhibit the results in the following three 
formulas: 


1 
| $2 (y)b*(y) dy = SmnklmAn" 
=4 


(33) 
+ cos u, COS u, (sech’ a + a’ tanh a). 
1 
| y o2(y)[a’o*(y) — o*’’(y)] dy = 2 COS Um COS Un(Am — An)? [DmAn (34) 
=~§ ¢ 
+ (20, — An)AwA, (2X, — 4° sech* a — 4a tanh a)], 
if m ¥ n. 
3 9) 
/ y'd*(y)[a"o%(y) — o%’’(y)] dy = d,(3-* — 27*u,”) 
7 (35) 


+ cos’ u,[\,2°'u, “a tanh a + A,u,* 
+ 8\7"(a’ sech’ a + a tanh a) — X,u,’a tanha — 4], 


Section II. Application to the question of stability and connections with the Orr- 
Sommerfeld equation. Truncating the system (17), we obtain. 


N 
tte = > key, m= 1,2,+++,N, (36) 


n=1 
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where N is a “large” finite integer and the constants k,,, are given as follows in terms 
of the previously introduced notation: 


Kenn = to(Cmn — bmn) — SmnXm¥- (37) 


Assuming that such an approximation is valid for sufficiently large values of N, it is 

clear that the stability of plane Poiseuille flow relative to perturbations of period 2ra~* 

with respect to x can be investigated by a numerical study of the roots o,(N), o2(N), 
, ay(N), of the determinantal equation, 


knn — Omae | = 0, mn=1,2,--- WN, (38) 


in «. In accordance with a familiar theory, if at least one of these roots has a positive 
real part, the truncated system is certainly unstable. A numerical study of these roots 
for different values of N also reveals the possibility of ordering them in such a way that 
the limit 
lim o,(N) = o, , [= 1,2,3,---, (39) 
ut 
exists for each fixed value of /. Evidently, if at least one of these limits appears to have 
a positive real part, one is justified in concluding that the Poiseuille flow is unstable 
under the given conditions, that is, the conditions under which a and v have fixed values. 
The numerical details will be described in Sect. ITI. 

It may be noticed that this discussion of stability avoids the Orr-Sommerfeld equa- 
tion on which previous discussions of stability have been based. That this avoidance is 
more apparent than actual may be seen in the following way. 

The well known Orr-Sommerfeld equation 


v(f'’’’ — 2a’ f’’ + a’ f) = ial(l — y® — o(f” — a’ f) + 2f] (40) 


is obtained by setting ¢(y, t) = e~***'f(y) in the partial differential equation (8) for ¢ 
and then multiplying by e**“* 

Evidently if, for some value of c, the Eq. (40) has a solution # 0 which satisfies the 
boundary conditions f(+ 1) = f’(+ 1) = 0, then we shall have instability if the imagi- 
nary part of ¢ is positive. The direct investigation of the values of ¢ for which this is 
possible has been the basis of previous studies on stability. This has been carried out 
by considering the general solution of (40), or approximations thereof, and then attempt- 
ing to determine the constants of integration so as to satisfy the boundary conditions. 
This leads to a homogeneous system of linear equations in the constants, the compati- 
bility of which system leads to a very complicated transcendental equation in c. 

An alternative way of treating this question is to expand the unknown solution of 
(40) in terms of our functions ¢,(y) previously introduced: 


fly) = > Pibi(Y), Dun = (fy On)- (41) 
It is possible to prove that such an expansion is possible if f satisfies the boundary con- 
ditions. Conversely, any function f(y) defined by (41) will satisfy the boundary conditions 
automatically (assuming convergence of the differentiated series), since each of the 
¢,’s satisfies the boundary conditions. Substituting the right member of (41) for f in 
(40), multiplying by ¢,,(y), and integrating with respect to y from — 1 to + 1, we find 








> 
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with the help of (11), integration by parts, and (18) that 
(—tat + An)Dmn = ia(—P + _ Cia) (42) 
n=1 


With the help of (37), this may be written in the form 


© 


Z. (Kmn + tach mn)/p, = 0, mo 2.2. +++, @, 


n=1 
‘ 


Truncation of this infinite system leads to the finite “approximate” system, 


. 
: (kmn + tacbmn)p, = O, 
n=1 

the compatibility condition for which is precisely Eq. (38) with o replaced by — ‘ac. 


The characteristic values c for the Orr-Sommerfeld equation may therefore be presumed 
to be identical with the o’s of (39) after multiplication by the factor ta™’ 

Section III. Discussion of numerical results. The majority of the numerical results 
obtained to date have been based on an eight-by-eight truncation of the system (17). 
Since there is nothing magical about this choice of order, a few words about the reason 
for it appear in order. Certainly on an a priori basis there is no reason to expect that 
such a low order system would be capable of representing the complicated behavior of 
the linearized Navier-Stokes equations quantitatively, for our method is essentially 
one of sampling, or a method of moments. On the other hand, if our method were to be 
of any help in understanding the mechanism of transition to turbulence, then it was 
desirable that a low order system contain the qualitative features of the problem of 
instability. Now, although we shall not do it here, it is easy to show that all two-by-two 
truncations of the system (17) are stable. In view of the alternate approach to this 
problem it is, moreover, reasonable to suppose that at least a fourth order system would 
be necessary in the method of moments before instability could be detected. The four- 
by-four truncated system (17) has been explored numerically on an 1103 Computer 
and no evidence of instability was found. Since we did not know when the phenomenon 
of instability would make its appearance, we decided to add equations and unknowns 
two at a time up to the limit of the code then available for complex eigenvalues, ten 
equations in ten unknowns, and hope for the best. We were fortunate in that this phe-. 
nomenon did occur in the six-by-six system.* 

At this time, then, it was decided that a complete curve of neutral stability for an 
eight-by-eight system would be computed, and that, in addition, curves of constant 
amplification would be of some help in indicating the regularities of the approximation. 
This has been done, and the results of these computations are compared with those of 
Lin [1] in Fig. 1. These results while most encouraging and of the right order of magnitude 
did not agree at all well quantitatively with those computed by Lin [1] and checked by 
Thomas [10]. After some preliminary results with a ten-by-ten truncation which did not 
indicate any major divergences from the eight-by-eight case, some preliminary evidence 
was found that the twenty-by-twenty system would give much better agreement with 
the results of Lin. The twenty-by-twenty truncated system of differential equations 
was integrated by use of the Runge-Kutta-Gill method and compared with a similar 


*The five-by-five system has not been explored. 
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integration of the eight-by-eight system. In the case of the eight-by-eight system, the 
integration was in agreement with the result of the eigenvalue calculations in that the 
system appeared either stable or unstable depending upon whether the corresponding 
eigenvalue has a real part which was negative or positive. While no finite run can be 
considered conclusive, a corresponding integration of the twenty-by-twenty system for 
a = 1.5 and» = 10“ indicated stability in agreement with the curve of Lin. This com- 
parison is indicated in Fig. 3. 

It turned out to be a major problem to develop machine codes and programs which 
would permit the calculation of the complex eigenvalues of the twenty-by-twenty un- 
symmetric complex matrices with sufficient accuracy. Fortunately this problem has 
recently been solved by the Computer Division of the Ramo-Wooldridge Corp. and it is 
possible to report a few numerical results at this time. The details and techniques 
necessary for these calculations will be reported elsewhere by members of the Computing 
Division. For the eight-by-eight case the computations were based on a routine pro- 
grammed at Convair, San Diego, Calif., and kindly furnished the Ramo-Wooldridge 
Corp. The routine was derived from a paper by Leppert [11]. The computations were 
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Fig. 3. Numerical integration of system of differential equations. 


TABLE 1 


Mazimum of the real parts of the eigenvalue 





N=8 

a\y 6 X 10-4 4x i 2x 10-* 10-+ 10-5 10-6 

0.7 —0.014652 —0.001594 —0.000159 
0.8 —0.003858 +0.007952 +0 .022386 +0.024191 
0.9 +0.009241 | +0.022902 +0.038748 | +0.040589 
1.0 —0.002491 +0.020052 +0.034649 +0 .051136 +0 .053007 
Ea +0.004910 +0 .028679 +0.043931 +0.060904 +0 .062810 
2 —0.008579 +0.010234 | +0.035119 +0.050919 +0.068348 +0 .070295 
1.3 —0.006618 +0.013419 +0 .039433 +0.055734 +0.073613 +0.075610 
1.4 —0.006915 +0.014613 +0.041783 | +0.058538 +0.076864 +0 .078909 
1.5 —0.009132 +0.014120 +0.042410 | +0.059539 +0.078274 +0.080370 
1.6 —0.012737 +0 .012323 +0.041588 +0.058961 +0 .078027 +0.080175 
1.7 —0.017162 +-0 .009597 +0.039580 +0 .057007 +0 .076278 +0.078473 
1.8 —0.021947 +0 .006241 +0 .036604 +0 .053832 +0.073121 +-0.075358 
1.9 —0.026806 +0 .002468 +0 .032828 +0.049534 +0 .068558 +0 .070828 
2.0 —0.001592 | +0.028380 +0 .044152 +0 .062457 +0 .064752 
3.1 | +0.023364 +0.037698 +0.054473 +0.056768 
2.2 +0.017894 +0.030233 +0 .043764 +0 .046007 
2.3 +0.012123 | +0.022052 | +0.027878 | +0.029702 
2.4 +0.006261 | +0.013936 +0.005404 | +0.000594 
2.5 +0.000555 | +0.006936 +0.001535 | +0.000155 
2.6 —0.004764 | +0.001538 | +0.000449 | +0.000045 
3:7 —0.002473 | —0.000118 —0.000012 
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carried out by assuming a given value of Reynolds number based on a channel half- 
width of unity and a maximum laminar velocity of unity, and exploring the system for 
different values of a. In this normalization the Reynolds number is numerically equal to 
the reciprocal of the viscosity. In none of the cases so far computed have we found more 
than one unstable eigenvalue for any given value of a, and our computations have 
furnished N distinct eigenvalues. 

According to Lin [1] the minimum critical Reynolds number as defined above is 
5300. The eight-by-eight system yielded an approximate minimum critical Reynolds 
number of around 2000, a result of the right order but quantitatively unsatisfactory. 
The currently available results shown in Table II for the twenty-by-twenty case imply 


TABLE 2 


Maximum of the real parts of the eigenvalue 


N = 20 

aw | 2x 10- | 15x10 | 10-4 

0.80 —0.000046 
0.85 | |  +0.002010 
0.90 —0.000774 +0 .003327 
0.95 —0.003076 +0.000540 |  +0.00377 
1.00 —0.002007 | +0.001012 +0 .003222 
1.05 +0.000540 +0.001562 
1.10 —0.000962 —0.001311 


that the minimum critical Reynolds number of this order of approximation is not less 
than 5000 nor greater than 6600. Linear interpolation of them indicates a value of 
approximately 5800. This value may well show even closer agreement with the results of 
Lin after further exploration of the nose of the neutral curve. It should also be observed 
that the minimum critical Reynolds number occurs in the neighborhood of a = 1, in 
general agreement with the results of Lin although the meager numerical results so far 
available make it impossible to decide the extent of the agreement. In view of the vast 
differences which exist between our techniques and those of the asymptotic theory the 
above agreement is remarkably good. 

Since the methods used in this paper are not asymptotic it is not reasonable to expect 
that the agreement given by the twenty-by-twenty approximation near the nose would 
hold uniformly for larger Reynolds numbers. In fact it does not. While nothing like a 
complete neutral curve has been computed in this case, isolated calculations for : 
Reynolds number of 10,000 show that the lower branch lies between a = 9.5 and a = 9 
while a similar calculation for a Reynolds number of 100,000 shows that the lower 
branch lies between a = 9 and @ = 8.5. Thus while this is evidence that the lower branch 
is monotonically decreasing it is also evidence that the lower branch lies above that of 
the curve computed by Lin. No corresponding calculations are yet available for the 
upper branch at large Reynolds numbers. This merely reflects the fact that our method 
involves series expansions which, while convergent for all Reynolds number, require the 
use of more terms the higher the Reynolds number. 

In spite of this inherent limitation we believe that we have furnished additional 
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evidence of the instability of plane Poiseuille flow and have indicated an alternative 
general method of attack on the problems of stability which is not only capable of 
complete mathematical justification but which involves only techniques familiar to 
almost all workers in the field. Last but not least is the fact that potentially these same 
techniques appear capable of furnishing a rigorous mathematical way of discussing the 
mechanism of transition to turbulence since in principle their scope is not limited to 
linear problems. Their success in the linear problem of this paper lends credence to the 
hope that it may eventually be possible to verify that transition to turbulence takes 
place as envisaged by E. Hopf [6]. 
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INEQUALITIES FOR EIGENVALUES OF SUPPORTED 
AND FREE PLATES* 


BY 
L. E. PAYNE 


University of Maryland 


1. Introduction. Much attention has been given in the literature to the classical 
vibration and buckling problems for a plate if the boundary is clamped (see for instance 
Weinstein [15], Aronszajn [1], Aronszajn and Donoghue [2], Polya and Szegé [11], 
Payne [9], and others). Many of the methods which have been proposed for handling 
such problems are not practical if the boundaries of the plate are supported, free, or 
satisfy mixed conditions. Of course one has a minimum principle for computing upper 
bounds for the eigenvalues, but in applications lower bounds are usually more important. 

In this report we establish certain “isoperimetric” inequalities for the eigenvalues 
in the classical vibration and buckling problems for supported and free plates. These 
inequalities involve the geometry of the plate and the eigenvalues of membranes having 
the same shape as that of the plate. 

It is well known that if the boundary of a supported plate is rectilinear the eigen- 
values in the vibration and buckling problems are related to those of the corresponding 
fixed membrane problem. If the boundary is not rectilinear this is no longer the case, 
and to my knowledge no relationship between the eigenvalues of the supported plate 
and those of the membrane is then known. We shall prove in this report that if the 
boundary is convex the eigenvalues of the supported plate are never greater than the 
corresponding eigenvalues of the fixed membrane of the same shape. 

Let us suppose that the plate occupies a region D bounded by a closed curve C in 
the xy plane. In comparing the eigenvalues of the plate problem to those of the membrane 
problem we assume the latter problem to be defined over the same region. We consider 


the following eigenvalue problems: 


: 9° 9” 
[. Au+ru=0 in D, u = 0 on (a = . s+ 2) 
OX oy 
: Ov 
lI. Av+p=0 in JD, —=(0 on C 
Ov 
II. A°w+ AAw =0 in D, w= 0, M(w) =0 on C 
IV. A*yg — Xe =0 in D, g = 0, M(e) = 0 on C 
¥. A°w + TAw =0 in D, M(w) = 0, Viw) = 0 on C 
VI. Aye —-y¢e=0 in D, M(ve) = 0, Vie) =0 on C 


where v denotes the outward normal on C and VM (W) and V(y) are defined as follows: 


*Received December 10, 1956. This research was supported in part by the U. 8S. Air Force under 
Contract No. AF 18(600)-573—monitored by the Office of Scientific Research, Air Research and 
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M(¥) = cAy + (1 — | 2 (2) ar, o (2%) au , O<e<} 


Ox/ Ov Ov \dy/ Ov 


" 0 0 0 fdw\ Ox 0 [dw)\ oy 
Viv) = — (Ay) — (1 — o) = : or a ee ae a 16 
Ov Os LOv \dy/ Ov Ov \ Ox/ Ov 
As is well known, I and II characterize the fixed and free membrane problems 
respectively. The buckling problem and vibration problem for a supported plate are 
given by III and IV, while V and VI are the corresponding problems for a free plate. 
In each of these problems it is well known that non-trivial solutions exist if and 
only if the constant A, u, A, 2, T and y take on certain discrete non-negative values. 


(1) 


We order these eigenvalues as 
a ae n> 'm, (2) 
the other eigenvalue being similarly ordered. We denote the eigenfunction corresponding 


to X,, by u, , ete. 
In this paper we establish the following isoperimetric inequalities for C convex: 


iE <= %, (3) 
os. Ay; (4) 
r. Ss (9) 
7% Sh (6) 
For arbitrary C we prove further that 

A, => (1 — o)pe + ad, , (7) 
\ fe es). (8) 

= 2 | 
A, > Tass > (1 — o)Miinsay21 » (9) 
m2 dA; , (10) 
Yar2 > Maly > (1 — open, ; (11) 
where the symbol | (n + 3)/2 |! is used to denote the largest integer contained in the 


quantity n+ 3)/2 


We first define 
- 0 0 p' i i ) 
X(y) = (1 — 0) p( 2%) + D(*)} +o | (Ay)? dA. (12) 
OX \ OY Ja 
D 
The eigenvalues in each of the plate problems which we consider are defined as the 
minimum value assumed by Wy) for a certain class of continuously differentiable 
functions ¥. In (12) D denotes the Dirichlet integral, i.e. 


Diy) = [| grad py GA. 


ve 


D 
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2. Buckling problem for the plate. A. Supported boundary. The eigenvalues A,, in 
problem III are defined by the following minimum principle: 


A, = min. A(y) (13) 


among all continuously differentiable functions which vanish on C and satisfy the 
relations 
r | av Ow OW ow , 
Wy, w,) I x < Pot hes « @ Pe Pee 
DW, w, di E Ox si Oy dy |. : , 2, sis 1, (14) 
D 


and 


Diy) = 1. (15) 


We note that if C is convex its radius of curvature p is everywhere positive. Likewise 
if we denote by r the distance from an arbitrary origin inside a convex curve C then the 
quantity h = r(dr/dv) is positive at every point on C. 

We assume first that the region is convex and choose for y, 


Y= Dau, , (16) 
s=1 


where the functions u,; are the eigenfunctions of the corresponding fixed membrane 
problem (I) and the a; are so determined as to satisfy (14) and (15). Insertion of (16) 
in (13) gives after an application of Green’s formula, 


= ‘ (1 — a) 0 = du; \’ = du, \’ 
i < A: Du) + = $ - ( dus) ( ous) 8. 7 
a da :D(u,) —" at da 5, + da. = ds. (17) 


We have used the fact that 
Diu; , u;) = 0, t x j. (18) 


Now the boundary integral in (17) may be rewritten as 


¢ : ( > a; out) a ( au)" ds = $ ¥ a; a ps a; Au; ds 
J ¢ i=] CO. i 
¢ 


a; 
by a4 Oy Ov fm 

¢. ou ] Tanne 
3 a ; (19) 

. ou;\ a 2 * Ox 

+$¢(da2 ) ° = —_ ds, 

J Gai Oy Os = . Ou; 

| a; pig 

. lions =| Op 


where 0/ds denotes differentiation along the curve C. The first boundary integral on 
the right vanishes if the region is convex since the u, satisfy 7. We now make use of the 
fact that on C 


OY _ Fg i 9E LS, Aus OY _ 
Ss dai Ox ast da, dy ds 0, (20) 


and find finally that 


F du; \? = du; \’ rites du;\ 
» > é s=-—p-|) ad . 9 
f Ov (: 1% Oy + ( = “ oy | . p p (> “6 “ap os 6. (21) 


c 
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Thus from (17) 


A, < >> aid; D(u;) < ry D> a?7Du) =», (22) 

i 1 4 1 
the last step following from the normalization (15). This establishes inequality (3). 
Note that we may obtain a better inequality for the first eigenvalue A, , for in this case 
we have only one eigenfunction u, in (16). The quantity on the right hand side of (21) 
can then be approximated. Let us suppose that the function is normalized according to 
(15) (then a, = 1, a; = 0,7 > 1). According to Rellich’s identity [12], 


a i y Ou; . aaa 
l =\X, [| u,aA = 1p n(%*) ds, (23) 


D " 
where h = r(dr/dv) and the origin is arbitrary. We choose the origin so as to minimize 
the quantity h,,,., on C and denote by A’ the resulting h,,,, . We obtain then the in- 
equality, 


9) — \ 
A, <x, - = 9 (24) 


As was previously mentioned the equality sign holds in (22) and (24) if the boundary 


C is rectilinear (pmax = ©). It is easily seen from (21) that this is the only case in which 


equality exists. 
Inequality (7) for : 
maximum-minimum principle for A, (see Courant-Hilbert [3]) or 


monotony principle (see Aronszajn and Donoghue [2]). One need merely observe that 


1 general plate (not necessarily convex) follows directly from the 
from the second 


in problem II, u,; = 0 (v, = constant) and x, is defined as 
; D(x) . 
He = min. —— = , (25) 
x dA 
/. 
for continuous functions x which satisfy the condition 
|| x dA = 0. (26) 
It follows then that, since the functions y in (13) vanish on C, 
€. wie) = 
[| dA I dh = 6. (27) 
JJ OF JJ oy 


D D 


Hence from (25 
ay ay)\’ (24) SI (24) a 
D{ Ox ) 2 mJ (24) =) : D oy = He oy 1A P (28) 


Inequality (7) then follows directly from the maximum-minimum principle. 
Inequality (8) is established by making use of the minimum principle for \, , i.e. 
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X, = min. — Dx) __ (29) 
| a 
D 


for continuous functions x which vanish on C and satisfy the condition 


[[ xu. =0, ee ee (30) 
D 


By Green’s theorem and Schwarz’s inequality we have 


_D&) FS» (31) 


= D(x) 
[| x dA ‘“ 
“s 
rf a a°x a°x ax \’ 
wore [23 (es) 6 
I] edliees J | oy” Ox Oy eA 


AGO 


Thus, 








A — P. = , Pe —_ 
5 D(x) D(x) +m . D(x) 
[/ (Ax)’ dA 
(x) (l—o)'p 
< 4. —& - — 29 
= Dix) 2 D(x) (32) 


lA 


,fi-e.f-«.... «f= « ... ee 
1 +( =*) +( 5) + +( 2)’ + | 


Thus we obtain for arbitrary x 
2 Hx 
. <2 MO) 


< T+ 6 DW) an 





We now choose x to be a linear combination of the first n eigenfunctions w, of problem 
III with the constants chosen to satisfy (30) and inequality (8) follows immediately. 
B. Free boundary. The eigenvalues I, of problem V are defined by Eqs. (13), (14), 
and (15). However, we now remove the restriction that y vanish on the boundary. It 
is clear that the first three eigenvalues will be zero. The corresponding eigenfunctions are 
w, = C,,w. = C, + Cyr, ws, = Cy + Cer + Coy, (34) 
where the C; are constants. 
We assume first that the plate is convex, and choose instead of (16) 
y= Didw,. (35) 
i=l 
Equation (17) follows with u, replaced by v,; as does (19). We have instead of (20), 
since dv,/dv = Don C, 


dy << Ov; Ox - Ov; Oy _ oes 
ov > a oe 2 “ay av = (36) 
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Instead of (21) we have in this case 
¢ 2 ( > Q; avs)’ + (x a a)" ds = - (> a; dus)’ ds < 0. (37) 
J dv — Ox) : oy 1 Os 
Thus from (17 
i > au, De Y<y, > a a; D(v;) = pnp « (38) 
<1 i=l 


This proves inequality (5). 

For the general plate (not necessarily convex) we shall now establish (9). The left 
hand side of this inequality follows from the minimum principle for I,,, mentioned 
previously in this section. (The eigenvalues I’, are defined by the right hand side of (13) 
under conditions (14) and (15). However, the trial functions y are not required to vanish 
on C). We note that any function which vanishes on C is orthogonal in Dirichlet norm 
to any harmonic function. (This is an immediate consequence of Green’s identity). 
Hence, if we choose for y functions which vanish on C, condition (14) is automatically 
satisfied for the first three eigenfunctions (34). In particular let 


y= > dw, (39) 


where the upper index s indicates that we are using the eigenvalues of the supported 
plate. The b; are chosen to satisfy (14) and (15). Then 


Ties S Uy) = z. bA(w;) = = b:A;D(w') < A, . (40) 
‘ 1 i 1 
In order to establish the right hand side of (9) we recall the definition for u, , ie. 
D(x) 
vw, = min. — x (41) 
I x dA 
i. 
for continuous functions x which satisfy the condition 
I] xsdA =0, i =1,2,+++,n—1. (42) 


We choose for x two functions 


da ow — . OW 
x = 8B ay X2 = > B; (43) 
> COL 2 


F Oy 


with the 6; chosen to satisfy condition (42). The w,; are the eigenfunctions of problem 


V. Then from (41 
D(x.) + Dow) + 72 |] | + aus | dA 
, L—od Ox oy 
Doxa) 4+. D(x2)_ us D 


D 





(44) 
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But since U(w, , w;) = D(w; , w;) = 0 for i # j this inequality reduces to 


a°| p(22:) + p(2”:) + ,%— I (Aw,)? aa | 
- L OX | oy l—o | 
n<d ——g——2——<-+-r.. (45) 
. . 3° D(w,) . 
We have then the inequality 
Psa 2 (1 7) Mn ’ (46) 
or equivalently 
re =~ (1 _ o) $2) 1043)/21 (47) 


which was the relationship to be established. 
3. Vibration problem for the plate. A. Supported boundary. The minimum principle 
for determining the eigenvalues 2? of problem IV is the following 


2; = min. A(y) (48) 


for continuously differentiable functions y which vanish on C and satisfy the conditions 


[[ wed =0, i=1,2,---,n-1, (49) 
and 
Hl y dA =1. (50) 
D 


Let us assume for the moment that C is convex. Proceeding exactly as we did in 
establishing (3) we prove with little difficulty that 


e<*. (51) 


Thus if we take Q, to be positive (4) follows immediately. One also obtains the improved 


9 an 
O< aa _ 2 — 2). (52) 


vy 


bound for Q , 


We return again to the case of a general curve C. Inequality (10) follows directly 
from the minimum principle for \, , (29) and (30). We choose for y 


y= Lew, (53) 
t=1 
where the c, are chosen to satisfy (30). (We use the g, of problem IV.) Then 


> ci; I g, dA 


9 ( 
1 MW Mg Fh gaa, (50 


7 D 
ff eas _ we Il ada 
D D 








118 L. E. PAYNE [Vol. XVI, No. 2 


which is (10). In a similar way we could have obtained the inequality, 


a a a O0O<m<n- 1. (55 


n 


From (55) and (8) we obtain at once 


(1 ) i 
- + <2. > See 0O<m<n-—1. (56) 


xX 2 
B. Free boundary. We turn now to problem VI. The minimum principle for y, is given 
by 
Y, = min. A(y) (57) 
among all continuously differentiable functions y which satisfy conditions (49) and 
(50). As in the case of the buckling problem for the free plate, the first three eigenvalues 
are zero with corresponding eigenfunctions 
ee k,, 
ig oe & Ee, (58) 
¢; =~ ky + ksx + key, k; = constant. 
For convex C the proof of inequality (6) follows along the same lines as the proof of (5). 
Hence we shall not repeat the details of the proof. 
For general C the first half of inequality (11) is proved in the same manner as was 


(10). The second half then follows from (9). To prove (11) we make use of the minimum 
principle for yu, given by (41) and (42). We may rewrite (41) as 


N(x) D(x) (59) 
X) x) 59) 


[f dA _— 


va 


D 


Hn => 


and choose for x 
x= 2D be, (60) 
i=] 

where the ¢, are the eigenvalues of problem VI. The 6, are chosen to satisfy (42), and 
in addition the conditions 

D(w, ,x) = D(w; , x) = 0, D(x) = 1, (61) 
where w, and w; are given by (34). [The condition D(w, , x) = 0 is satisfied identically.] 
Then 


* by; I/ w, dA 
bie ‘D__ - 


Mn S (62) 


lA 
2 
zw 
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which is the first half of (11). In a similar way we establish the inequality 

Ya+3 = Slecudd ome ’ 0 < m < n— - (63) 
With the help of (49) this inequality gives 


Yn 3 => (1 x 7) Mo + mbt n—m+3)/2| » 0 < m < » 1. (64) 


Tr 


The right hand side is a maximum if m = n — 1. We obtain then 
Yass > (1 a 7) Moby +1 ° (65) 


4. Concluding remarks. The eigenvalues of the vibration and buckling problems 
for a supported or free plate have been bounded in terms of the eigenvalues of the 
corresponding fixed or free membrane problems. If these cannot be determined ex- 
plicitly they are still more easily bounded than the eigenvalues of the plate problems 
which we have considered. For instance the eigenvalues of the fixed membrane are 
monotone with region while those in the plate problems which we have considered are 
not. Under certain symmetry conditions one can obtain a lower bound for the first 
non-zero eigenvalue in the free membrane problem which involves only the geometry 
of the region (see Payne and Weinberger [10]). We have also the isoperimetric inequalities 


> r' , GY 2.4048, 
. (66) 
wn <p 1.8412, 


where A is the area of D. The first inequality was proved independently by G. Faber 
[4], and E. Krahn [7], and the second was formulated by Kornhauser and Stakgold [6], 
proved under certain restrictive conditions by G. Szegé [13] and established for a general 


region by Weinberger [14]. From (66), (8) and (56) we have the inequalities 
_(l+o)" (l+o) |[?xj - 
A, 2 _ _ , Q, > Ez _ (67) 


We note further that inequalities (3) to (10) hold in any number of dimensions. 
The right hand side of (9) must merely be changed to read 


Pyea-1)+2 2 (1 — ob, (68) 


where N is the number of dimensions. 
In N dimensions inequality (11) is replaced by 


Yatn > Maly se > (1 os 7) Moby . (69) 


In case of a vibrating free plate in the shape of a square, Nakata and Fujita [8] 
have applied the method of Kato [5] and after considerable computational work have 
obtained both upper and lower bounds for the first non-zero eigenvalue. For ¢ = 0.225 
and a square of side 2 they obtained the bounds 


3.418 < ¥, < 3.554. (70) 


However, their method is not completely rigorous. If we make use of the symmetry of 
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the region’ we obtain the bounds 
3.06 < y, < 4.94. (71) 


The upper bound is, of course, not good and a close bound can easily be obtained anyway 
by the usual Rayleigh-Ritz procedure. However, one finds immediately, with no com- 
putational work whatsoever, the lower bound given in (71) which in application may 
be quite good enough. 
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A FREE BOUNDARY VALUE PROBLEM FOR THE 
HEAT EQUATION* 
BY 
W. L. MIRANKER 


Beil Telephone Laboratories, Inc. 


1. Introduction. In this paper we will prove the existence of a solution of a free 
boundary value problem for the heat equation. We will accomplish this by demonstrat- 
ing the existence of a solution to a non-linear integro-differential equation. 

Let D be the domain 0 < t,0 < x < R(t), R(O) = A, indicated in Fig. 1. The boundary 


value problem is 








t 
D x= R(t) 
—_- x 
0 oe 
Fic. 1. 

u, = a’u,,(2, HeD, 
u.(0, t) = —G(t) Gi) > 0, 
u[R(t), t] = T. = constant, (1) 


u(z,0)= F(a) >T O< 2x A, 
u,{R(t), t] = B — CR,t) B>O,C > O are constants. 


Differentiation is denoted by a subscript whether it is a partial derivative of a function 
of two variables or an ordinary derivative of a function of one variable. 

This problem with A = 0 has been discussed by several authors, see for example 
(1, 2, 6], but thus far no existence proof has been established**. The problem describes 
the physical phenomena of evaporation, fusion, sublimation, etc. For example with 
B = 0, (1) could refer to the following situation. A long metal rod insulated at the sides 
has begun to melt at one end (x = 0). The layer of liquid metal is A units deep and 
has some initial temperature distribution, F(x). The critical temperature 7’. is the melt- 
ing point of the metal. At z = 0 heat is applied to the rod at a known rate proportional 


*Received March 26, 1957. 
**For our proof it is essential that A ~ 0. See, however, the remark at the end of the existence 


proof below. 
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to G(t). As the process continues the interface, R(t), between liquid and solid advances 
down the rod. 

This problem is called a free boundary problem since the part R(t) of the boundary 
of D is unknown. The additional boundary condition (5) which would over-determine 
the problem were R(t) known compensates for the free boundary. 

We set B = 0 to achieve clarity in our discussion, and we introduce the transfor- 


mation 


a’ =¢e-—T,, 

R’ = (A*/a)R, (2) 
g =G, 

f =F, 

D’ = D. 


Upon dropping the primes the free boundary value problem then becomes 


U, = Uz, , (2, HeD, (3) 
u(0, )) = —g() < 0, (4) 

u[R(t), t] = 0, (5) 
u(x, 0) = f(x) > 0, (6) 
u,(R(i), t] = —R Ad. (7) 


We require that g(t) and f(x) be continuous and have the following additional proper- 
ties: g(t) is differentiable; f(x) is continuously differentiable for 0 < 2 < A; f(x) = 
f(— x); f(z) = Oforx > A;f,(A) < 0; and f(0) = — g(0). An f(x) with these properties 


18 


f(x) = me x — g0)|x!+ g0OA — A’," lz| <A, 
f(x) = 0, otherwise. 

2. Method of solution. Our method of solution will be to apply the method of I. 
Kolodner [4] and derive a functional equation for R(t). We will show that the existence 
of a solution with certain properties of the functional equation implies the existence 
of a solution to the free-boundary problem. We will solve the functional equation by 
the method of contracting maps (Picard iterations). 


The method of Kolodner: Let x = p(t) be a continuously differentiable function and 
such that p(0) = A. Consider the function u’(z, t) defined as 
w= +w’ + f’, (2.1) 
where 


v= —ir”” [ (t — 7) °'*p,(r) exp {—[3la — pla) ](t — 1)” 7} dr, (2.2) 
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we = —he [ (t — 1)"2p,(r) exp {—[Fla + o(Dt — 1 "P} dr (2.3) 


+n? | (t — 7) '°g(r) exp {—[3a(t — 7) ']} dr, 


and 
t 


f? = (rl) '” | f(g) exp {-[}(z — Bt? P} dr. (2.4) 

This function is a solution of the heat equation in the domain of Fig. 1 with p(t) 
replacing R(t). We will now calculate the same boundary and initial conditions of u*(z, t) 
from 2.1 which are prescribed for u(2, ¢) in the statement of the problem (3)-(7). That 
is by computing uv’ from 2.1 and then letting (2, t) approach the boundaries of the domain 
of Fig. 1 we shall obtain the boundary conditions in question for uv’ and u? . 

The arguments used in doing this are lengthy and are based upon well-known formulas 
for the integral solutions of the heat equation (see [3]). 

The results are 


ue(p(t) + 0, 2) — ub[p() — 0, t] = pf), (2.5) 
ui(0+,t) = —g(d, (2.6) 
w(x, O+) = f’(x, 0+) = f(z). (2.7) 


(2.6) makes use of the evenness of f(z). 
(2.6) and (2.7) show that u’(z, t) satisfies the boundary and initial condition (4) 
and (6) which are required of u(x, t). (2.5) almost does the same for (7). If in (2.5) we 


require that 


u(p(t) + 0, t) = 0 (2.8) 
we will have 
ue[p(t) — 0, t] = —p,(2) (2.9) 


which is the condition (7) for u. We will see that (2.8) is an integro-differential equation 
for p(t). To show that a solution to this integro-differential equation furnishes a solution 
to the free boundary problem, we need only apply Green’s formula, 
2 I| (u,)? dx dt = $ u’ dx + Quu, dt (2.10) 
D 


to the domain D, of Fig. 2 and to the function w’. 
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Using (2.7) and the fact that f(z) = 0 for x > A and also using (2.8) we have 


2 [| (uz)? dz dt = 2 | u*(a, thut(a, t) dt 


f 


J p(t) 


(2.11) 


(w)"(e, 4) de — | “(w’)"[p(t), t]p, dt. 


If we now let a — ~ and observe that (2.1)—(2.4) imply that u’(a, ¢) and u2(a, t) > 0, 


then (2.11) becomes 


2 Hl (u2)* dx dit + | (uw’)"(x, tj) dx = — | (u’)*[p(t), t]p,(t) at. (2.12) 


Deo 


If we require that 
> 0 (2.13) 


p 


then the right-hand side of (2.12) is not positive and (2.12) shows that vu? = O for x > p(t) 
I \ 


and all ¢. This in turn implies that u’(p + 0, t) = 0, and since w’ is continuous along 
x = p(t) that u’(p, t) = 0. 

Thus the requirement (2.13) makes yu’ satisfy the condition (5) for wu. We see that 
if p(t) is a smooth function for which p(0) = A and p, > O and such that u2(p(t) + 0, 0) 


= 0, then the function uw’ solves the boundary value problem (3)—(7). The condition 
(2.8) yields from (2.1)—(2.4) the following integro-differential equation for p 


p,(t) = pl(i)r '”” [ airi(t — Ne exp {—[Fp(A(t — T ldr 


— 4," l p-(r)(t — 7) **[p(t) — p(r)] exp {—[3[e() — o(d](t — 7)? )} dr 


” (2.14) 
— p(r)(t — 7)**[p(t) + p(7)] exp {—[3[o() + o(z)](t — 1) P} ar 
+ 4n '?¢*" [ f(é)[p(t) — &] exp {—[4[p() — Et} | dr, 
p(0) = A. 
We introduce the abbreviation 
p(t) = F(p, p. , 9, f, A, ) = F(p) (2.15) 


for (2.14). 
If we can solve this equation and if its solution p(t), has a non-negative derivative 
then p(t) is R(t), the free boundary, and w’ is a solution of the boundary value problem. 
3. Properties of »7(t). In this section we deduce some properties of the integro- 
differential equation (2.14) and of the free boundary which will be of use in our existence 
proof. 
Lemma 1. If r(¢) is a continuously differentiable function for ¢ > 0 then 


p(0) = lim F(r) = —f,(A). (3.1) 
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Proof. The proof proceeds according to the methods of [3]. The first three integrals 
on the right in (2.14) tend to zero when ¢ tends to zero while the fourth tends to — f,(A). 

Lemma 2. If R,(é) exists and is continuous then R,(t) > 0. 

Proof. From (6) and (7) we see that R,(0) = — f,(A) > 0. Suppose to the contrary 
that R,(t) > 0. Let t’ be the smallest positive value of ¢ for which R,(t) = 0. Let D,- be 
that part of D where t < t’. u, is a solution of the heat equation in D,. . Since g(0), f.(z), 
and R,(t) are continuous and f,(0) = — g(0) and — f,(A) = R,(0), u, is continuous on 
that part of the boundary of D,- which is not on the line ¢ = ¢’. By the maximum princi- 
ple [5], the maximum of u, occurs on this part of the boundary of D,. . Since — g < 0, 
f'(x) < 0, and R,(t) < 0 for 0 < t < ?#’, this maximum must occur at z = R(t’), t = t’. 
Thus u, < 0 everywhere in D,. . Thus since u = 0 for z = R(t), we have that u > 0 in 
D,. . Then in the closure of D,, , every point on the free boundary is a minimum point 
of u. Now at a minimum point of u, the outward drawn derivative in a characteristic 
direction must be strictly negative*. That is, u, must be strictly negative along the free 
boundary. But at ¢ = t’, u. = — R, = 0. This contradiction implies the result. 

4. Existence. To demonstrate the existence of a solution to the free boundary 
problem (3)-(7), we must show that the integro-differential equation (2.14) or (2.15) 
possesses a solution p(t) for which p,(t) > 0. To do this we will use the principle of con- 
tracting mappings. We will first obtain the existence of a solution p(t) with p,(t) > 0 
in the small and then using Lemma 2 of Sec. 3, show that this solution exists for all 
t>0O 

Existence in the small. Let B be the Banach space of continuously differentiable 
functions {p(t)}, 0 < t < T for some fixed T to be specified. Let the norm be 


p(t) || = | p(O) | + lub | p,(d |. (4.1) 
O<stsT 
Thus convergence in B is uniform convergence of the function and its continuous de- 
rivative. 
Let G be the closed set in B whose elements satisfy the following properties 
(a) p(0) = A, 
(b) p.(0) = — f(A), (4.2) 
(c) O<l<p()<K, 


where 1 < — f,(A) and K > — f,(A) are to be specified. 

Now consider the map p? = F(p). We will show that if p, and p, ¢ G, then p{ and 
pheGand |! pf — pi || < C(T) |! p: — po ||, where O < C(T) < 1 for T sufficiently small. 
These statements show that p/ = F(p) is a map of G into G which is continuous and 
moreover contracting. Thus by the principle of contracting mappings pj = F(p) has a 
fixed point in G. This fixed point is our solution and possesses by [4.2(c)] a positive de- 
rivative. We now proceed with the proof. 


If p(0) = A and p,(0) = — f,(A) then the same is true for p’(0) and p/(0). The first 
since p} = F(p) is an integro-differential equation and we may arbitrarily require 
o/(0) = A and the second is the assertion of Lemma 1 of Sec. 3. 

In Appendix 1, we show that 

Il er — p21] < C(T) || op: — 2 ||, rr s poe G, (4.3) 


*This is an unpublished result due to L. Nirenberg. 
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where C(7) = C(T, 1, K, A, g, f) < 1 for T sufficiently small. In Appendix 2 we show 
that p/ is continuous. In Appendix 3 we show that 


| p(t) — [—f.(A)] | < Ce”, (4.4) 


uniformly for pe G. 
Thus we fix 7 so that C(7’) < 1 and 


C,T’” < max {| l — [-—f,(A)]|,| K — [-—f,(A)] ]}.- 


Our map is then into and contracting and a solution p(t) of (2.14) exists up to time 7’. 
Moreover this solution has a positive derivative. 

Existence in the large. The proof of existence in the large proceeds as the above proof. 
The Banach space is now B, the set of continuously differentiable functions p(t), 
0 <t< 7, where 7, > T is to be specified. We use the same norm as in B. G is replaced 
by G, , those functions in B, which up to time 7’ are equal to the solution R(¢) which we 
have just shown to exist. A is replaced by R(T) and — f,(A) by R,(T). l and K are 
replaced by two other numbers J, and K, withO < 1, < R,(T) and K, > R,(T). 

With this setup the requirements of the principle of contracting mappings are 
satisfied here in essentially the same way as above. Thus we produce a 7, such that 


for 7, — T > 0 and sufficiently small, we have existence of R(T) up to time 7, and 
moreover R,(t) > 0 for 0 < t < T, . We see that we may iterate this procedure and 
produce a sequence of 7’; ,2 = 1,2, --- , such that R() exists and R,(t) > Ofort < T,. 


There remains only to show that 7; — © . From the form of the estimates in the appen- 
dices and the definitions of the sets G; we see that this will be the case if R,(t) never 
vanishes for then we may always extend our solution slightly further. But the vanishing 
of R,(t) is ruled out by Lemma 2 of Sec. 3. Thus 7; — © and F(t) exists for all time. 

Remark. We have mentioned that for our proof it is essential that 4 + 0. A passage 
to the limit as A — 0 is suggested to obtain the solution with A = 0. This limit procedure 
would be legitimate if for some sequence R*(t) with A tending to zero, the slopes, R?(t), 
have a positive lower bound. For in this event the set of functions R*(t) are equi-con- 
tinuously differentiable and a simple compactness argument justifies the passage to the 
limit. If the condition (7) is changed to u, [R(t), t] = B — R,(t), B > O, then an obvious 
extension of Lemma 2 yields a positive lower bound for the slopes R?(t) and in this 
case our method yields existence for the case A = 0. 

Acknowledgment. The author is grateful to Professor I. Kolodner and to Professor 
J. Keller for their assistance in preparing this paper. 


Appendix 1. Continuity and contracting of the map pi = F(p). In this appendix we 
show that the map p/ = F(p) is continuous for p e G is contracting. Let u(t) and v(t) 
be two generic elements in G. Let ui = F(u) and v} = F(v). We have 

u’ — v’ || = max |! F(u) — Fv) || 
t 7 
since u’(0) = v’(0) = A. 


The computation of the right-hand side of (A1) is divided into four parts corre- 
sponding to each of the four integrals occurring in F. The computations for the first 
three integrals are so alike that we illustrate the computation involving the third and 
fourth integrals only. 

(i) THtrD INTEGRAL. We are led to consider 














NI 


1958] FREE BOUNDARY VALUE PROBLEM FOR HEAT EQUATION 12 


at 


I,+I1,= 390°” | u(r)(t — 7)~**[u(t) + u(r) — v(t) — v(7)] 


-exp {—[}fu(d) + u(r) (t — 7) '7} dr | 


4 Ag? | [ (v.(7) — u(r) }o(t) + o(s)\(t — 


-exp {—[3lo(d) + o(7)](t — 1)" F} dr) 
By the law of the mean, we have 
I, < 4Knr'” [ | u(t) + u(r) — v(t) — v(r) | (t — 1)? | eo — 22’) | dr, 
where z is [3[p(t) + p(r)] (t — r)~'”*] for some p « B. Thus 


I< Kier? \lu—vl] [ @-— 9? |e" — 2) | ar. 


it 
“0 


Now 
2A + It - AD tot) A+KE 5 
at — 7 = 9 4 SG “O . 
Then 
L< Kir? \lu oll | © 2 + A + KOE = 7] 
-exp {—[}(2A + l(t — 7)7'?)} dr. 
Let 
g 1(2A + l(t — 7)”, 2do = 3(2A + lt)(t — 1°? dr. 
Then 
I, < LAr y= F | “(2A a lt) ' | ee [l a 8o°7(A + Kit(2A + lt)~?] de. 
J2A+1t/2e'/2 
For 7, we have 
I, <2 u—v | (v(t) + v(r))(t — 7)~*? exp {—[3fo() + v(t — 7) "7 PF} dr 
S* u—v|| (A + Kb) | (t — 7) *? exp {—[3(24 + l(t — A? PY dr 
» At Kt .* 
tn ~Nieere I ant de 


(4 1/2 : - _— 
< i(—) u—v||\(A + KO(2A + I) exp {-32(2A + Wt}. 
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ii. FouRTH INTEGRAL. For this integral we have 


nA | 
I1+1,=340°°t" | f(£)[u(t) — v(t)] exp {—4[u(t) — Et") dé 
e A ' 
| nA 
+4! ff ( — ellesp {Huo — eFC 
| /—A 
— exp {—3[u(t) — é]'t"}] dé | 
nA 
I, < 3)” ||[u-0| | | {© | exp {—3[u(t) — E]Pe-"}? dé 
4 
nd A 
Lint)? ll u — | + [ | #(€) | exp {—2[u(é) — E]’t"}? dé 
J—A a 
I,-+1,. 
where \ = Xi 
IT, < 3(wt)"’? ||}u —v || (A — 2X) max | f) 
A<t<h 
I; < 3(rt)"'’*(A — d)* max | f,(€) | (Al) 
A<t<x 
since f(.1) = 0. 
I, < 3rt)'?(\ + A) exp {-12(A + lt — d)*t"} max | f(é |. (A2) 


—ASESK 
From (A1) and (A2) we see that if as ¢ > 0 X(t) > A faster than t'”* but slower than 
1/2 
t’’~ then 
cs — i s uo | const o(1) 

The computation involving J, is similar to the one just conducted and will be omitted. 

Appendix 2. Continuous differentiability of an image under F, In this appendix we 
show that if p? = F(p), p « G then p? is a continuous function of ¢t. The computations 
are slight variations of the computations in Appendix 1. Therefore we will carry them 
out only for the first integral in F(p). 

Let a > 0 and consider 


a 


rola + Ab) | g(r)(a + At — 7)” exp {—[3e(@ + AD(a + At — 7) I dr 
—T pla | a(r)(a T “exp {—[p(a g¢— “"P tl dr 


= 7? o(a + Af) | a(r)(a + At — 7 exp {—[}p(a + At)\(a + At — 7)7'”*]*} dr 


+” | g(r) [pla + Ad)(a + At — 7)” exp {—[3p(a + Ad(a + At — 7) J} 


— p(a)(a — 7)” exp {—[3p(a)(a — 7)" }] dr 
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In what follows, we use the fact that p e G implies that p(a), p.(a), p(a + Ad), 
p.(a + At) > 0. 
Consider first J, . Let ¢ = p(a + At)/2(a + At — r)'”. Then 


I, sie 4y_7 3? [ ge~*" dx, r= 4 o(a + At)(At)"’”? 


vz 


< 4Ma-”? [ ede, M = max | g(é) |. 


O<t<T 


Since 
i mgt 1 -z? 
ec“ ds<sz-e', x>0 
* 2x 


we have that 
I, < 4Mn7'”(At)'”?/p(a + Ad). 


For J, we have on application of the law of the mean of the differential calculus: 


I< ram | {| exo ~—<—o + 3 p(i(t — 7%" 


0 \ 


\-1/2) 1 -1/2 1 \-3/2 
— p(i(t — 7) {3 p.(t)(t — 7) + 4 p(t)(t — 7) 


“exp | -[3 p(t(t — ol} dr, 


where the quantity in the curly brackets is to be evaluated at some value of ¢ in the 
open interval 0 < a < t < a+ At. The integrand is thus finite and the integral exists. 
Appendix 3. The initial value of p{ . In this appendix we show that | p/ (i) — 
[— f,(A)] | < ©€,(0)'”. We sketch the ideas of the proof since the computations are 
variations of those in Appendix 1. 
We note first that the first three integrals vanish as t — 0. For the first integral, 


a”? o(7) [ g(7)(t — 7)*” exp {—[2e((t — 7)" P} dr, 


the singularity at £ = 7 in the exponential [p(0) = A > 0] causes the integrand to be 
bounded in the range of integration. This bound depends on the minimum of p(t) eG 
and this is bounded below by A + IT. Thus the entire integral vanishes with its upper 
limit uniformly in G. 

For the second integral, 


te? [ p.()folt) — p(t — 7)? exp {—[Blo() — ala) t — 7)7F} dr, 


the term [p(t) — p(r)]/(¢ — 7) approaches p,(0) = — f,(A) uniformly in G. Thus the 
integral tends to zero like t'’”* uniformly in G. 
In the third integral, 


er o.(r)[p(t) + p(r)](t — 7) °? exp {—[3[o(d) — p(x) ](t — 7)7'”7]"} dr, 
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the singularity in the exponent causes the integral to vanish with its upper limit as in 
the first integral. 

We have, referring to [3], observed that the fourth integral tends to — f,(A) ast 0. 
By examining the proof of this process one may observe that the limit is approached 
like t'’’ as t > 0. 
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AIR DRAG EFFECT ON A SATELLITE ORBIT 
DESCRIBED BY DIFFERENCE EQUATIONS IN THE REVOLUTION NUMBER* 


BY 
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Abstract. For a satellite on an orbit which is affected by air drag, difference equa- 
tions are derived whose solutions express changes in orbital size and shape and give 
the satellite’s time behavior as functions of revolution number. The results are obtained 
from a first order perturbation theory using a small air drag parameter. 

Introduction. The behavior of a satellite orbit under the effects of air drag have 
been discussed by Petersen [1]', the present author [2], and others. Using the method 
of variation of parameters, I derived some approximate expressions for the decay of 
eccentricity with radius, the decay of radius with revolution number, and the growth 
of time with revolution number. Those results appear to be especially useful for eccen- 
tricities which are quite small, say « < 0.01, although a complete investigation of their 
range of validity has not been investigated. 

An alternative approach is available which may be preferable for larger eccentricities. 
It consists of a perturbation method in a small parameter, carried only to the first 
approximation, followed by the derivation of an exact relationship between values at 
the beginning and the end of any orbital revolution. Thus one obtains a set of difference 
equations in the values of the problem's dependent variables which occur at each full 
orbital revolution. 

Because the approach leads to difference equations, it appears well suited to numerical 
treatment. It offers a significant advantage even when a high speed digital computer is 
available, for it obviates the numerical integration around each orbital revolution 
which can impose a significant burden on machine storage capacity and result in a 
prohibitively long solution time. It is clear that in a conventional numerical solution 
in which the satellite is followed step-wise around each revolution, even the finest 
practicable integration interval and the greatest care in handling errors may not avoid 
a significant accumulation of error over tens of thousands of revolutions. On the other 
hand, the difference equation method developed here allows the errors for any complete 
revolution to be made as small as desired without excessive difficulty. 

Outline of the method. Like [2], this treatment begins with the drag law and equa- 
tions of motion discussed by Petersen [1]. The following notation is used: 

r radial distance from center of earth to satellite 

8 angular advance of satellite from arbitrary interially-fixed reference line in 

plane of orbit 

K product of earth mass and constant of gravitation 


*Received March 28, 1957; revised manuscript received August 27, 1957. 

1Numbers in square brackets refer to the bibliography at the end of the paper. 

2These apply only to the case where the aerodynamic force acts only as a drag along the line of the 
velocity vector. Lift forces are assumed absent. Other assumptions are given with the notation below. 
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Cp drag coefficient, assumed constant 

satellite projected area on plane normal to velocity vector, assumed constant® 
m_ satellite mass 

a equatorial radius of earth 

h_ altitude, r — a (neglecting oblateness of earth) 

p(h) standard stable atmospheric density function 


In these terms, the equations of motion are 


r= 1B me BE ative? + 8, (1) 
B+ ors = Ses p(h)r8'(r'? + 1°78"), (2) 
It is convenient to introduce auxiliary dimensionless variables 
E=a Fe (3) 
7 = Ka (r*B’)’, (4) 


and to regard them as functions of 8. Primes henceforth denote differentiation with 


respect to 8 as the new independent variable. Also define the parameter. 
y = CpAa/m. (5) 


The equations which are basic in the sequal then become 


Y+E= 2, (6) 


(agé — Q@) eee = 

n n nun . (1 + (& ¢)”)' ° (4) 
é 

A key assumption is made, as in [2], that the density function is locally exponential. 

That is, in a restricted neighborhood of any particular é-value, the function p(h) can be 


approximated by 


pak —a) =e - (8) 


where, in general, a and 6 vary slowly with the £& (altitude) level at which p 


is evaluated. 

The following procedure is used. 

1.£(8) and 7(8) of Eqs. 6 and 7 are expressed as power series in v (second and higher 
order terms of which are subsequently neglected). 

2. The zeroth approximation, elliptic motion, is found. 

3. The first approximations to &(27), £’(27) and n(27) are found in terms of &(0). 
£’(0), and (0) as series in increasing powers of eccentricity, under the single assumption 


that p(£) can be approximated by an exponential function over the range of £ repre- 


sented in one revolution, as described above. 
4. These values of £(2), &’(2r), n(2r) are used as initial conditions for the next 
revolution, and the same procedure is followed. 


’This is the case if the satellite is spherical, and is closely true for other vehicles whose attitude 


relative to the earth is kept nearly constant by a suitable control system. 
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5. The results are precisely the same for connecting the initial conditions for the 
1)st revolution to the initial conditions for the jth revolution. 
6. The difference equations are returned to the original physical variables of interest. 
Perturbation solution. Let ¢,, ¢’, n; (j = 1, 2,...) represent values of £, d§/d8 and 
n at the beginning of the jth revolution and é;(8) and 7;(8) represent the functions 
£(8) and (8) during the jth revolution. Let &; = a/r, , & = — aVy/r,Vy and m = 
Ka/r{Vz be given initial conditions, in terms of initial radius, horizontal velocity, and 
vertical velocity. However, there is no loss in generality in supposing that 8 is measured 
from a point at which £{ = 0. 
Assume that £ and 7 can be expanded in power series in a small parameter (vp), 
where p, is any value of density such that p(h)/po is of the order of unity in the altitude 
region of interest. In particular, it may be very convenient to choose pp = p(at;* — a). 


This implies the representation 


£,(B) = £,°(B) + (vp E(B) + (vpo)ES?(B) + -: , (9) 
n(8) = ;°'(3) + (vpo)n}''(B) + (vpo)?n)?(B) + °°: . (10) 
The zeroth approximation during the jth revolution satisfies 
e778) + £8) = n\"(8), (11) 
n, (8) = 0. (12) 
The first approximation satisfies ; 
&°°'"(B) + &:"(8) = 95"(8), (13) 
0 = Nt + [Eo lt as 
As initial conditions, take 
£, = £52 (Qr) + (vp)é,"\(2n) + -:- 
Ef = E22) + (mpad (2x) + +++ Gj = 2,3, +) (15) 
n, = ni 2(2m) + (vpo)n (Qn) +--+ J 
and 
(0 E, + (vp,)0 + ---, £70) = 0, 7,10) = » + p)O+ ---. (16) 
We begin with the first revolution, solving Eqs. 11 and 12 to obtain 
m (8)=m, (17) 
£, (8) = (1 + «, cos 8), (18) 
where 
e, = (&/m) — 1. (19) 
Using these results in Eq. 14, 
n'(8) = m f pin.(1 + €, cos 8)] Be coed a dg. (20) 
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Now the assumption discussed in connection with Eq. 8 enables us to simplify the 

the integrand of Eq. 20 by 
p[m(1 + e, cos B)] = p(n.) exp (maye, Cos B). (21) 


The remainder of the integrand can be expanded easily in a power series in ¢, the resulting 
coefficients involving the Legendre polynomials. These, in turn, can be developed in 
Fourier cosine series, whence 


l+ 2cosBt+e)" _ 5 T | (22 
errer , = 2d ( e) T',(€) cos np, (22) 
where 
ZZ 3: 4 
m=1+ e+ = 
‘eae Tat + 
r=1+ oY 9 3 4 m (23) 
1 8 € a 39 * 
1 l 9 2, 
2 r + i6f + 


An explicit general expression for T, can be given, but only in a rather complicated 
form. For n > 3, the constant term (independent of e) in I, is 


] : (Qn — 2k — 3) "kk + 1) 
> ies a ‘(n — k)'"(n — k — 2)! 


Combining the results of Eqs. 21 and 22, Eq. 20 becomes 


f -) af 
p(7:) = ; ‘ 
m" (8) = — > (—e,) T,(e,) | exp (7,@,€, cos 8) cos np dB. (24) 
The integral in Eq. 24 does not have to be evaluated because we need only n,'' (27), 
and by [3], p. 181, 
[ exp (x cos 8) cos n8 dB = 2rI,(x), (25) 


I,,(x) being the usual notation for the Bessel function of the first kind for purely imaginary 
argument. Thus we have 
_ 2mp(m) 


m (2s) = nt > (—e,)"T,,(e,)/,,(niaye,). (26) 
Po n=0 


Next we must solve Eq. 13. The solution can be written 


af 
1 ° 1 7 
&; (B) = | sin (8 — 2x)n, (zx) dx (27) 
which, in view of Eq. 24, becomes 
ded < 8 “a 
¢,'’(8) = “ > (—e)'T,.(e) | sin (8 — 2) i exp (n,a,€, cos y) cos ny dy dx, 
) Qn 0 J 


nB 


(mn) < 
a > (—e,) T,,(e;) | exp (a, 7,€, COs y) cos ny [cos (B — y) — 1] dy 
Da Geo 0 
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From this it follows that 


2 (7) , , / 
&, (2) = a" 7 2, (—«,) T,.(€:) (Zne i (mares) + (n/moae, — 1)Z,.(mai€,) ]. (29) 


Similarly, 


me p\ 7) an 
é y= — P » (-«) r.(é,) 


When 8 = 27, the integrand of Eq. 30 becomes an odd function of y with period 2. 
Therefore, £;'’ (27) = 0 and one concludes that, to first order in »v, the second cycle 
begins with the same kind of initial conditions as the first cycle, i.e., conditions appropriate 
to apogee or perigee. By induction, £ vanishes identically in j. 

Equations 26 and 29 relate conditions at the end of the first cycle to those at the 
beginning of that cycle. Moreover, because £/ is zero for all j, exactly the same argument 
can be used for any cycle provided the eccentricity and a-value appropriate to that 


8 
| exp (,a,€, Cos y) cos ny sin (8 — y) dy. (30) 


70 


cycle are used. 
\ set of interim results complete through first order terms in the perturbation param- 


eter Is 
t &. + 2rvp(n;) ) x (—e,)"T.(e,) [L,4,:(n;aje;) + (n/n,ja;e; — II,(nja;e;)], (31) 


' 


20 


Misr = 0; t+ 2rvp(n;) » (—e;)"T(e;)7,.(n,a@;€;), (32) 


n=0 


re 


, 2rvp(n;) , a : . 
e, + ; 2. (—e,)"T,(e,) [T,..(n,aje;) + (n/njaje; — 2 — €;)1,(njaje;)]. (33) 
Equation 33 follows from the fact that ¢,; is defined as &;/n; — 1 for all j. Also, for j = 1 
we have the initial conditions given at the beginning of this section. 

Some remarks on accuracy. The accuracy with which one can find £, 7 and ¢ from 
Eqs. 31-33 depends principally upon the accuracy with which v and p are known, and is 
therefore a physical problem. However, these equations do involve power series in e, 
so it is of some interest to see how the errors arising from truncating these series are 
related to the number of terms required. 

In analyzing the orbit of artificial satellites, it should be noted that e« < 0.2 is probably 
almost certain, even ¢ = 0.1 being rather large. Also, one might typically be interested 
in altitudes to within one mile (¢ about one part in 4000) at the end of 20,000 revolutions. 

The latter is suggested by some of the numerical work in Ref. [2].) Because errors in 
the right-hand terms of Eq 31 are additive, this means that for each 7 the term must 
be correct to about one part in a million. With « < 0.2, noting that all of the leading 
coefficients of the T,, are of the order of unity, it follows that one generally must carry 
e’ terms in the I, . This accuracy normally will be adequate in Eqs. 32 and 33 as well. 

Not all the T, terms, of course, are equally important in contributing to the error, 
for the coefficients of T’, in the infinite series (which converge) eventually become small. 
Unfortunately, it is not easy to establish just how fast these series converge, where 
they may be truncated, or how the accuracy requirements on the I’, may be decreased 
with increasing n. However, some indication can be obtained by noting that the product 
nae may be about 25 for « = 0.2. (A value of about 12 for « = 0.1 is given in [2]). The 
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asymptotic expansions of J, and J,,, can be used in Eq. 31, and it is found that about 
twenty-five terms of the series are required before the residual drops below 10™°. 

This argument must be done separately for each equation in view of the accuracy 
requirements one may wish to impose on the basic dependent variables. I believe, 
though, that this last result may be considered typical. 

Thus the number of terms required to give high accuracy, in effect, to make the 
method adequate for a much greater eccentricity than that permitted by earlier methods, 
is substantial. However, it generally will be considerably easier to carry this number of 
terms in a numerical computation than to do a piecewise numerical integration separately 
around each orbital revolution, and to continue this through 20,000 or more revolutions. 

Radial and time behavior. It is desired to return to the original physical variables 
of interest, r and ¢t. From r = a/¢ it follows that to within terms of the order of v’, 


) 


or 


“a@M™V oa, vy site " . ‘ ‘ 
ae = r; p(n;) DR (—e;)"T.(€;) [Tnsi(njaje;) + (n/a;e; — 1)1,(9;0;€;)]. (34) 
Also 
dB (K\” #(8) 
# _ (E)" $0. 33) 
dt a, 7 (8) 
Integrating dt/d8 over one revolution as was done in finding 7;'’(6), 
4 1 \3\1 of a2" 1)/\ 2,1 > 
ri(1 + e,) 2r mm, (8) 2m, (8) ) 
a. (nil ) a oe | Wh VL _ Sts: V2 | de>. (36 
2 1 \ K 1 a é)° 2 - Jo Qe) 2/3) £ *(B) J ) 


Each of these integrals can be evaluated, albeit tediously. Doing this, and generalizing 
the argument to the jth revolution, 


3 \ 1/2 (= — 
P 7 Vv 3) € 
t.,=tc+ 7 is] E t s i " (Qr) 
i "\K(1 — €;)" = n; | 2 bs (37) 
of 
e(1+¢é)" < pas ; ) 
+ ae ZL = r’.(e;) (7, 1(9,;@;€;) + Re nace)I} | 


Equations 31, 32, 33, 36 and 37 form a complete set of difference equations from 
which the change of radius, shape and time with revolution number can be inferred. The 
n; enter this set only as auxiliary variables, and are not concerned in the final geometry 


of the orbit. 
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1. Introduction. Properties of Gaussian processes must be investigated in many 
problems involving the analysis of random noise. Almost as common as the normal 
distribution is the Rayleigh distribution which occurs in work on radar, the detection 
of signals in the presence of noise, properties of a sine wave plus noise, ete. [1, 2, 3, 4]. 
It is the purpose of this paper to investigate certain properties of generalized Rayleigh 
processes. 

We shall use the notation X/a, ¥(7)] to indicate that ¥ is a Gaussian process with 
mean a and autocorrelation function (7). [That is, ¥(7) = ({X(d) — a] [X(t + 7) — a), 
where Y(t) is a member of X.] A Rayleigh process, 3t, may then be defined as 


W = Xla, Wr) + Mlb, W(>)], (1.1) 


where ¥ and ¥) are independent Gaussian processes. Certain classical results that are 
associated with 9{ may be found in the above quoted references. For example, 
(a) The first order probability distribution of M: 


| PWR) = = exp [—(R® + A’) 2yolta( 24) , (1.2) 
where y ¥(0), A* = a’ + Bb’ and J, is the Bessel function of the first kind and order 


zero With purely imaginary argument. 
b) The second order distribution of 9: 


| RR, R? + R? | ( ARR: ) 
R, ,R.) = => exp | — =e A SS) > 1.3) 
P ve Fi —- exp | Q¥(1 — *) J Nod — 2°) “a 
where \ = y(r)/y is the normalized autocorrelation function and we have assumed 
a 0 = b. 
(c) The correlation function of 9: 
C(r) = y[2E() — (1 — d*)K()], (1.4) 


where AK(\) and E(A) are the complete elliptic integrals of the first and second kind 
respectively. 
In the present paper we wish to consider generalized Rayleigh processes which we 


define as follows: Let X,[a, , ¥(7)], X¥:[a. , W(r)], --- , ¥xvlaw , ¥(r)] be N independent 


Gaussian processes. Then we define the generalized Rayleigh process, vt, as 
+2 +2 + 271/2 
R= (Ki + Xs + oes + Ay]. (1.5) 





*Received April 8, 1957; revised manuscript received August 2, 1957. 
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The general problem is to compute the joint M/-dimensional distribution of 3. We 
have not been so fortunate. While it is of course possible to write down the distribution 
as a multiple integral, our objective has been to obtain compact usable formulas. Towards 
this end we have calculated: 
(i) The first order distribution of 9: 
N/2 > 
p(k) = £ (“) exp [—(R? + A’) /2ypo]/ r-nn(E4) : (1.6) 


where A’? = a? + a3 + --- + a and J,(z) is the Bessel function of the first kind and 


order v with purely imaginary argument. 
(ii) The joint probability density of 9: 


re (N-2)/2 
p(R, , R,) = 2RsHe TNA) (2) exp {—[(RE + Rn + Au} 





—— 3 
Yo N— 2 wv (1.7) 
foe T/9) bo 
x be an( _" : 2) (2MaR.R,)IiloR,)1(wR,), 
k=(N/2)—1 4 a / 
where 
; —— (1.8) 


"Tew? ** ai -d) 


and (°) is the binomial coefficient. 
(iii) The correlation function of 9: 








riwe2r/ana/ NW + 1) 
2y,(1 po d”) N naya(N + 1) ie 7 : 
Pic 27 (Mtl, N41 X y) ay 


I'"(N/2) 


where 2/', (a, 8, y, 2) is the hypergeometric function. 
(iv) The three dimensional distribution of ®: 


am BBE ME ay | — Mul + Mel + Malt) 








R, ,R,,R;) = - at pe : 
‘ | (M,.M,,M;,)""~”” 2|M| 
a (1.10) 
/ k+1—(N/2 k N gr k ~~ 2 la 
x p> : (—1) os ( i yd A ) H@hahG ; 
where 
RMyR,  ,_ RMaR, | _ RMuk, } 
a — M | ’ B _— M | ? oS | M | ’ (1.11) 
|| ¥) Wr) Wr +8) | 


M=||\%) vO vw | (1.12) 


Wirt) v® vO 


is the correlation matrix, | M | the determinant of M and M,, the cofactors of M. We 
have assumed A = 0, that is, the unbiased case. 
Other auxiliary equations, relations and special results will be pointed out as we 


develop the above formulas. Note that the results of (a), (b), and (c) are all subsumed 





tw 
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under (i), (ii), and (iii) above. There are, of course, many other ways in which one 
could generalize Rayleigh processes. Some of the methods and results are described in 
the references [5, 6, 7,11]. 

The authors wish to thank the referee for pointing out certain techniques which 
resulted in simplifications in the proofs of some results and led them to generalizations 
of other theorems. 


2. One dimensional biased Rayleigh. Our first problem is to find the fr.f. of R, 


R = [(Xi+ X;+--- + Xx)”, (2.1) 
where XY, , X, , --- , Xy are independent random variables with the biased Gaussian 
distribution 

] al P 
ae so, ee “aera - ; 
p(X, (Qnv,)? exp [—(X, a;)/2pol, 1 Pe 3 , N (2.2) 


Introduce an orthogonal set of coordinates in N-space with the vector A = 
fa, , @, , «++ , @y} as the polar axis. If @ is the angle between R and A, then the com- 
ponents of R in a spherical coordinate system in N-space are 


R, = R cos 6 


a-2 
R. =Rsin 6 [|sing, cos¢a-}, @ =2,3,---,N-—1 (2.3) 


k=1 


N-2 
Ry = Rsin 6 [[ sing . 


k=1 


Unit vectors in the directions XY, , X, , --- , Xy form an orthogonal system in N-space 


which are related to R, , R,, --- , Ry by an orthogonal transformation. Thus in particular 
N N 

R?= > R2 = > xX? (2.4) 
a=l a=l 


and the Jacobian from R, ; R, yee Ry to R, 6, eg * hoy y-o is identical with the 
Jacobian from X, , X,, --- , Xy to R, 6, & . Note that ¢y_, ranges from 0 to 27 while 
6, and ¢, , --: , éy-3 range from 0 to =. 

The fr.f. of R is given by the marginal distribution 


p(R) = | doyaf ++ [ Jp(X., +++ Xv) dO dd, dbs +++ doen, (2.5) 


where the Jacobian ./ is 


N-3 


J = R*' sin? 6 [| sin”~*~* ¢, . (2.6) 
a=l1 
Since the X, are independent, 
N 1 N 
p(X,,°*:,Xy) = I] px) = — =o GS | - > (X; -— a)’ 2. (2.7) 
i=1 (2ryo) r i=l 


But 
> (X; — a)? = > (X? — 2a,X; + a) = R’? — 2A-R + A’, (2.8) 
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where A -R is the inner product of the vectors A and RF and 
A-R = AR cos 0. (2.9) 
Hence 
p(R) = 2nR — fe WIL | [ sin*~? @ exp (+ AR cos 0/Y) dé 
(2m po) = (2.10) 
. 


x [] [ sin’ *~* ¢, do, 


Now the first integral is 


(1\_(N — 1\(2W%\"~?” AR’ 
r(5)r(*5 (24) Tex-ayal =) 
and the product is 
Tp a+ 1 | T : 
0 a(S 2°’ 3) ~ PUN — 13/2) 


Using these results in Eq. (2.10) yields Eq. (1.6). 
It is easily seen that as A approaches zero, that is, as a, , --- , @y approach zero, 
p(R) as given by Eq. (1.6) approaches 
' 2R* 
p(R) = —————mmnm—= @XD (—"/2y¥p) (2.11) 
(2yo)"’ “T'(N/2) 


which is the square root of the chi-squared distribution. 


Note that if NV is an odd integer, say N = 2n + 3, then Eq. (1.6) may be written 
in terms of elementary functions by use of the identity 
: (2a i se ad” sinh br 
a 1/2)() -_ 17/2 : 2) | z 2.12) 
ry d(x") x 


3. The two dimensional distribution of 9. We shall establish a formula slightly 
more general than that given by Eq. (1.7). The proof is no longer nor more difficult 
than a direct derivation of this equation. 

Let 

R,=([(Xi+Xst+e+ Xut " R, = [Y: + Veter: + Vu)”, 
where the joint 2N-dimensional distribution of X, , --- , Xv, Y: , «++ , Yw is of the 
form [[*_, p (X, , Y,) and 
p(X; , Y;) = [2 (det M)'’*]" exp {— 4[C,,(X; — a,)° 
+ 2C,(X; — ai)(¥; — b) + Cr( Ys — b,)"}}. 


The covariance matrix is 


Wy M,, M,.| 7 | Var X, Cov (X, , Y.) | 
| Mas Me | | Cov (Y;, X;) Var Y, | 
and C = || C;, || = M™~’ is the inverse matrix. 
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We shall show that 


p(k, , R,) ORR FR? Yk) 9° VOL, RDLG@R,) LARC wR, & , (3.1) 


k=(N/2)-1 


where 


2-9 “T((N — 2]/2) _ 


(det M)*”?(w,C ow») _— 


? 


S = C,,(R? + | A |?) + C.(R? + | BP) + 2C,.4-B, 


@) | CuwA + C..B |, eo. = | CB + Cid 1) 


6 Cee) «1 (cos a) = Gegenbauer function. 


The quantities A and B are vectors with components fa, , a2, +--+ , dy} and {b, , b., 
- , by} respectively while bars denote the magnitude of a vector: for example, 


N 
| CuA + C.2.B |? = DO (Cra; + Cid)’. 


We have let a be the angle between the vectors C,,A + C,.B and C..B + C,,A. Also, 
C*(t), the Gegenbauer function, is the coefficient of 8" in the expansion of (1 —2 Bt + 8”)~’ 
in ascending powers of 8. 


If we let A = B, My, = Mo. = Yo, Mio = Ay , then Eq. (3.1) reduces to Eq. (1.7). 
Note that in this case a = 0 and 


a T(2v + n) 
" = . 
ol) = ~ (Qn) 





If, further, we let A = 0, then p(R, , R,) assumes the elegant closed form 


pi, ,R RR 
(5 vera. — x) 
-exD | - _RE+R |r (_KR, (3.2) 
PLT Qo = SO? WoT = 9) 


One application of Eq. (3.2) is to a generalization of the second example in the paper 
by Barrett and Lampard [8]. This was pointed out to the authors in a private communi- 


cation from Dr. Lampard. 
We now turn to the evaluation of the density function of Eq. (3.1). To obtain 


p(R, , R,) we must evaluate the integral 


- NV 
AR, ,R,)=[ — a¥ [ dX T[ p(X. , Y)) 
? |} Ry J|Xi\=Re i=) (3.3) 


= [2r (det M)'?]*e"*” [, exp [(Ci2A + C2.B)-¥] d¥ 


Ry 


| exp [(Ci,4 + Ci2B — Ci2Y)-X] eX, 
J |X\|=R, 
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where X and Y are the N-dimensional vectors with components {X, , X., --- , Xv} 
and |Y,, Y.,-°-: , Yy} respectively. The notation f,x;~», means: integrate over the 
surface of an N-dimensional sphere of radius R, . 

To evaluate the integral with respect to X we choose the direction of Z = C,,A + 
C,.B — Ci2Y as the polar axis, make an orthogonal transformation and then a general 
spherical coordinate transformation as in Sect. 2. We find that 





| ee" dX = (2aR)" | ZOO Tw-nall Z | Re), (3.4) 
where | Z |? = w? + C?,R? — 2Ci,(CiA + C2B)-Y. Thus p(R, , R,) becomes the 
multiple integral 
p(R, » Be) = [QrR-' det M] N 2,-S 2 

(3.5) 


[exp (CA + CoB) ¥) | ZO? Lyn nll Z| R,) AY. 
JIVi=Ry 

Now the vectors C,,A + C.,.B,C,,A + C).B and Y span a three dimensional subspace 

of N-dimensional space. Choose the direction of C,,A + C,,B as the z-axis and let the 
zx-plane be the plane spanned by C,,A + C,.B and C,,A + C,..B. Let a be the angle 
between these two vectors, @ the angle between Y and C,,A + C,.B, y the angle between 
Y and C,,A + C,.B and ¢ the angle between the projection of Y on the xy-plane and the 


x-axis. Then 
cos YW = cosa cos 6 + sinasin 6 cos ¢. (3.6) 
Making an orthogonal transformation and then a general spherical coordinate trans- 


formation as before we have 


p(k. ,R,) = — <2 = Dr * ee —) Pei | exp (w,R, sin a sin 6 cos ¢) sin‘ “¢ do 
mEN — 2}/2) \2 det J Je 


x | exp.(w,R, cos a cos 0) | Z er oe v~-9)a(| Z R,) sin* “9dé0 


cy rai Me eee (3.7) 
tae (R ) ‘( R, ¢,)" = en 
3 \27/ \det M W> SIN a/ 
4 | exp (wR, cos a cos 6)I;y-3)/2(wokt, sin asin 6) sin'*~"”79 
| ees ee ie oe a OR 
where | Z |* = wi + Ci,R? — 2C,.w,R, cos 6. Using the generalized Neumann addition 
theorem [10] 
r —(N/2) +1 wip 9 V-2)/27 'N 2 2 : k 
(R, | Z I.y-2)/2(| Z| R.) = 2 ni) 2 KD 
\ _ k N/2 1 (3.8) 


Gea" Lee Rye“ Lee, RCS aa (0 6 


and the formula [9] 
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| exp (y cos a cos 6)I,_1,2(y sin a sin 6)C%, (cos 6) sin’*’”?6@ de 
J0 a 


on 1/2 
é (**) sin’~”"aC (cos a)I,..(y) 


the formula of Eq. (3.7) reduces to Eq. (3.1). 
4. The correlation function of Jt. By definition, the correlation function C(r) of 8 


is given by 


ci) =| | RRpR, ,R,) dR. aR, , (4.1) 


where the density function p(R, , R,) is given by Eq. (3.2) and we recall that 
p(R. , R,) = Ofor R, , R, < 0. If we perform the integration, 


Oy, (1 — ,2\(Nt2072 2 2p r 
a. ICS FE. (n+ v 4 1) 








IN ms of sf 

r(3) ‘pine +3) 
7 (4.2) 
an) 

eaieaeen 2 (N+1 N+1 NN ., 
= Qy(1 — “YP? — ey (Xt? N41 Ny), 
2 

where 2F,, is the hypergeometric function. 

If N is even, say N = 2», then by use of well known identities [9] involving the 


hypergeometric function 


_— \2\"+1 v—1 Ay 4 "| 
Qy(1 — r*)"** a’ |e K(a) (4.3) 


C(r) = : wa —s < ~~, 
- re) day 


a-xPr w1-] 

In particular, if VN = 2, » = 1 and C(r) reduces to the formula of Eq. (1.4). Since the 
derivatives of the elliptic integrals K(A) and E(A) are linear combinations of K and £ 
with rational functions of \ as coefficients, we note that C(r), for N even, can also be 
expressed as a linear combination of K(A) and E(A) with rational functions of A as 
coefficients. 


If N is odd, say N = 2» + 1, then 
) = a rs oe API, hs ; , .’). 
, a(n | 

r(5)r( +5) 


fy 1 ) ra d — \2\-1/2 ‘iis __ \2)-1/2 
2FA(I, 1, 5, ¥*) =f 10 - xy aretan A — XY"? 


The identity 


implies that in this case C(r) is an elementary function of X. 

5. Three dimensional unbiased Rayleigh. We shall derive Eq. (1.10) for a more 
general covariance matrix than that given by Eq. (1.12). 

Let 


R; = (Xj, + Xin + +++ + Xin)’, =f = 1,2,3 











144 K. 8. MILLER, R. I. BERNSTEIN AND L. E. BLUMENSON [Vol. XVI, No. 2 


where the joint 3N-dimensional distribution of the X,; is of the form tk 1 WX; ; 
Ast, As;) 20 

W(X1;, Xo; , X3:) = [(2n* (det M)]"'” exp (—3X/CX,). (5.1) 
The covariance matrix is M = || M,; || and C = || C;; || = M~* isits inverse. The column 
vector X; has components {X,; , Xo; , X3;} and X{ is its transpose. 
We shall prove by induction on N that 


px(R, ,R.,R;) = 2| M|~*”(R,R.R;) exp E 7 | 
] I 
Kx Yd vuhi@l(aly), (5.2) 


where 





and a, 8, y are defined in Eq. (1.11). 

The proof of Eq. (5.2) for N = 2 follows the pattern used in Sec. 3 and will be omitted. 
Of course one could also use the methods of that section to obtain Eq. (5.2) directly 
and thus entirely avoid the induction proof. However, the induction proof, in some 
respects, is neater and will perhaps suggest additional generalizations. 

Let us therefore assume the validity of Eq. (5.2) for N. We shall prove it true for 
N + 1. Perhaps the easiest way of doing this is by the method of moment generating 
functions. Let u; = R? , 7 = 1, 2,3 in px(R, , R, , R;). Then 


' EE. & ibis ast 
Du(U; , Ue, Us) = = ~ exp | -3 be C, «| ‘> v,1,(a’)T,(B’) 1,y’), (5.3) 
- i k N/2)-1 


where a’ = uj’7C,,u}’”, etc. Using the fact that the integral of py(u, , U2 , v3) over the 
whole range is equal to one, that the C;; appear only in the exponent in Eq. (5.3) in a 
linear fashion, and setting C;; = C/,; — 2t; we find, after integrating and dropping the 


primes that the moment generating function of py(u, , U2 , Us) is 
twit,&,)@(Cr"*\|c-arTr, (5.4) 
where 
été OF O 
T=ti0 & © 
Oo @ fz 


Equation (5.4) is valid if t; < 4C;,; . 
Now let 


S, R: a x; ver = Uz; + x; wee 9 j = 1.2.3 














1958} GENERALIZED RAYLEIGH PROCESSES 145 
where py(u, , U2 , Us) is given by Eq. (5.3) and p(Xi.wai » X2.wer » Xs,we1) is given by 
Eq. (5.1). In order to prove the induction we must show that the moment generating 


function of S, , S. , S; is equal to 
Mysilt , te ’ ts) = my(t, , te ; ts) | C pre | C — 27 lo 


Let $(t; , tg , fs) be the moment generating function of S, , S, , S; . Then 
O(t, , te , be) = malt, ,& , &) [ 


~ o 3 
| | exp | > Xe} 
J —@ — oo —-@7 i=1 


‘D(X war ’ Xo,n+1 ’ X3,w+1) dX, vei AX» v+1 AX 3 noi 


/2 


ult, 4 WIC" |c -~ar, 
which proves the induction hypothesis. An identification of 17 with the matrix of Eq. 
(1.12) establishes Eq. (1.10). 
Note the similarity in form of the first order biased distribution, Eq. (1.6), and 
the second order unbiased, Eq. (3.2), as well as the similarity of the second order biased, 
Eq. (1.7) and the third order unbiased, Eq. (1.10). 
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BOOK REVIEWS 


An introduction to mathematical economics. By D. W. Bushaw and R. W. Clower. Richard 
D. Irwin, Inc., Homewood, IIl., 1957. xii + 345 pp. $7.00. 


As stated in the preface, the purpose of the book is “to provide an introductory account of mathe- 
matical methods of economic analysis which is accessible to persons with a limited training in 
mathematics.” Part I (‘‘The Economics’’), contains a very readable account of the theory of price 
determination; its scope is indicated by the following chapter headings: Introduction—Macroeconomic 
Statics—Macroeconomic Dynamics I: Isolated Markets—Macroeconomic Dynamics II: Multiple 
Markets—Microeconomics I: The Theory of Consumer Behavior—Microeconomics II: The Theory 
of Business Behavior—Concluding Generalizations. Part II (‘“The Mathematics’’) contains material 
on functions, graphs, and equations; the elements of calculus; quadratic forms and determinants; 
difference equations and differential equations; and maxima and minima. While the book aims primarily 
at introducing economists to the mathematics that is useful in their field, the first part would serve 
equally well to introduce applied mathematicians to the field of mathematical economics. 


W. PRAGER 


The analysis of structures. By Nicholas John Hoff. John Wiley & Sons, Inc., New York 
and Chapman & Hall, Ltd., London, 1956. x + 493 pp. $9.50. 


Theory of structures is presented in this book mainly from the point of view of the principle of 
virtual displacements and the theorems of minimum potential and minimum complementary energy. 
Linear and non-linear elastic materials are considered and there is an introduction to plastic limit analy- 
sis and to inelastic buckling. Various topics treated on conventional lines are included, e.g. moment- 
distribution and slope-deflection methods for continuous beams and frames. Mathematical topics are 
developed from scratch; nothing beyond elementary calculus is assumed, and sometimes not even that. 
More knowledge of structures is occasionally presumed; e.g. results are quoted without proof from the 
theories of thin plates, of St. Venant torsion, and of non-uniform torsion. 

The book is in four parts, the first giving an exposition of the principle of virtual work and its direct 
y\blems. The second part presents the theorem of minimum potential energy 
ised to derive differential equations and as the basis for methods of succes- 


application to structural pr 
and shows how this can be 

sive approximation (Rayleigh-Ritz and Galerkin). Part Three deals with buckling problems, using the 
theorem of minimum potential energy to derive differential equations of adjacent equilibrium states, 
tayleigh-Timoshenko procedure. Part Four presents the theorem of minimum 


and also illustrating the 
and includes an introduction to plastic limit analysis. Throughout the book 


complementary energy 
the structural problems chosen for illustration are taken from both the aeronautical and civil engineering 
fields. 

The book should be of enormous interest and usefulness to those working in advanced structural 
analysis. The main defect of the book, in the reviewer’s opinion, is the author’s tendency to add dis- 
cussions of subsidiary topics to such an extent that the main physical and mathematical ideas do not 
stand out clearly and strongly. As one example, the author’s treatment of plastic limit analysis methods 
is based on rather full discussions of several examples of limit load calculations, using collapse mech- 
nisms and virtual work. But, although the ideas are all exemplified, statements of the two fundamental 
theorems (furnishing upper and lower bounds, respectively) are not explicitly given. The reviewer 
believes that the whole book could have been made more readable and more useful as an introductory 
text by sharper separation of basic concepts from illustrative and explanatory material. Still, these are 
matters of style and pedagogy; what is most impressive about the book is the obvious mastery of the 
author of both the practical and theoretical aspects of the subject, and his delight in it. 


P. S. Symonps 


(Continued on p. 154) 


147 


HEAT FLOW IN A CYLINDER* 
BY 
WALTER P. REID 


Michigan State University 


Summary. A formal, analytical solution is given to the problem of finding the tem- 
perature in a cylinder which is gaining or losing heat by radiation exchange with a 
concentric, thin-walled metal tube. The tube in turn is radiating to, or receiving radiation 
from, its surroundings. A general treatment of systems with this type of boundary 
condition is presented first, and then the given problem is solved as a particular case. 

Introduction. The problem being considered is that of radial heat conduction in an 
infinite, right circular cylinder surrounded by an air space, and then by a thin-walled, 
concentric, cylindrical tube. The tube is assumed to be a good thermal conductor, so 
that its temperature may be considered to be uniform throughout at all times. The rates 
of heat exchange between cylinder and tube, and between tube and surroundings are 
taken to be proportional to the respective temperature differences. The temperature 
of the surroundings is assumed to be some specified function of the time. 

The feature that makes this problem of interest is the boundary condition. It is 
given by Eqs. (31) and (32), although the nature of the boundary condition is more 
easily seen if w (the temperature of the surrounding metal tube) is eliminated from 
those equations. The surface temperature of the cylinder, u, is then seen to satisfy the 
equation 
o,u + a2 ~ + oa, f + oa, a = ¢;9(t), (t+ 0), 
where the o’s are constants. This same boundary condition was handled earlier [1] for 
a semi-infinite medium. However, in the present case, the medium is of finite extent in 
the direction of the heat flow, and so the method of Ref. [1] cannot be used here. 

It will be seen that the answer to the present problem is obtained as a series of func- 
tions which do not satisfy the usual orthogonality equation, but which satisfy instead 
the modified orthogonality relation given in Eq. (14). Before tackling the given heat 
flow problem, it is convenient to derive this modified orthogonality condition and some 
other preliminary results. 

Modified orthogonality relations. Let 


Pe, seth, ont, ( 


£ | Pe ou | + (QM) + ESM]u, = 0, (2) 
or or 


4) : 
v, = (c, + D, 2 VL: es (3) 





*Received July 8, 1957. The solution to this problem was begun while the author was at the U. S. 
Naval Ordnance Test Station, China Lake, Calif. 





148 WALTER P. REID [Vol. XVI, No. 2 


/ 4 ; . 
W, (Cc, + D, = WA » U), (4) 
, 


Then 


fa) Ou ; Du, a) a) a) i, Oe = 
2 | Pon, Mt — Pom | =u 3 [PM] —w2[rm Me] — @ 
Or or or or or or or 
(£ —_ £7) S(r)u,u 2 (8) 
So 
(A. — Ce (. — (' £ & — £7) | NS(r TP dr 
f 9u(L. ,t Ju, (1 t 
= (A, — C,82)(A. — C82) P(L,)| u(L, , ) 2 - a (L, , 0 
. r or or 
. ; oul L, . t) Metis . % 
— P(L,)| u,(L, , t Nx ) an (L, . t Me | (9) 
or or 
ra i i t) ) 
«fs, ~ CDP) —— a CH+a,-D ran { 
O7 7 or 
oul(L, , t) a 
— (A. — C.&)P(L ———- Eg — C,&.) + (B, — Dé) : Ju a (10) 
: or or 
yi / t) > \ } 
= (A, — C,#)P(1 (4 + B, - ) — a(c, + 2) Ju t) 
O7 or ¥ ” OF 
Ou (I t 4 , ) 
— (A, — Crt?) P(L,) “SS | (4 +B, 2) — 2(c, + D,*) jut ) (1) 
6] or] or 
Ou Pay - 
wit, — Cher.) = L (2 — £)1 
or 
( 4 t) 
~ (A, — CaDP(L,) 2 gt (12) 
or 
Hence 


e-2) 1, — CEA — Cat) | SQ) ar 


Sie du,(L, , t) 
+ (A, - CAP), 


or 


Au,(Ls , t) 
— (A, - C,#2)P(L,)w, — — ay = 0. (13) 


or 


The expression in brackets is thus equal to zero for n ¥ k. This gives a modified ortho- 
gonality relation. In the special case when C, = C, = D, = D, = 0, the development 
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above is seen to reduce to the usual orthogonality treatment, as given, for example, 
in Ref. [2]. Equation (13) is not useful if A, and C, are both zero, or if A, and C, are 
both zero. However, other equations may be derived in a similar manner for these cases. 
For example, by multiplying Eq. (8) through by (B, — D,&) (A, — C.&) and then 
proceeding as in Eqs. (9)—(13), one obtains: 

Lo 


(& — &, [ B, — D,&)(Az — Cok) | S(r)u,u, dr 


du,(Le 9 t) 
0 


— (B, — D,t)P(L.)w, 


— (A, — C.82)P(L, vu(Ly , 0 | =0. (14) 


This equation will be used in solving the present problem. 
An expression for the /‘th coefficient. If 
= R,(r)T,(?), Ur = dnl’, (t), w, = ¥,T',(d), (15) 


) and (15) one finds 


T 0| us — Dg)(A — Cat) | ~ S()R? dr — P(L)¥,(B, — D,#) Ue 
Jr, or 


then from Eqs. (14 


aL 


~ iat A. cL) | = (B, — D,t?) |  S()uR, dr 
JL, 


2, 0R,(L.) it _ ' ' 
— P(L,)u(B, — Di&) —o P(L,)v( A, — C.é)R,AL,). (16) 
Both sides of this equation were obtained by summing the brackets of Eq. (14) over n. 
On the left side of Eq. (16), however, use has been made of Eq. (15), and the fact that 
the brackets of Eq. (14) is zero for n ¥ k, while the right side of Eq. (16) was determined 
by interchanging operations of integration and summation, and using Eqs. (1) and 


15). 
If one sets 0 in Eq. (16), he obtains 7,(0), and hence the kth coefficient of the 


Fourier series for uv. For the case when the surroundings are at constant temperatures, 
Iq. (16) therefore enables one to complete the solution to the problem. Equation (16) 
with ¢ = 0 is usually obtained by getting the Fourier series for the initial value of u, 
multiplying by the kth eigenfunction, and integrating. The orthogonality relation then 
causes all terms of the series but the kth to drop out, and thus gives the kth coefficient. 
However, that procedure is not convenient to use in the present case because of the 
nature of the orthogonality relation—see Eq. (14)—hence the need for Eq. (16). 

The fact that the temperature of the surroundings in the present problem of heat 
conduction in a eylinder is a prescribed function of the time can be taken care of by 
means of Duhamel’s theorem. However, this turns out to be a bit lengthy in the present 
case. For this reason, the following development will be used instead. 

A substitute for Duhamel’s theorem. Let it be assumed that, in addition to Eqs. 
the following equations are also satisfied 


} 


(1)-(6) and (15), 





a , ou S(r) du ™ 
or | Po) me + u(r) - x ot? (17) 
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, i » fi d v(t) . 
(4, — Cie) + (B, — Die) 2 ula 51 ,t) = £,G,(0 - (4 + ct) — (18) 
Pa re 2, O 7 d o\ w(t 
| (4s — C.£,) + (B, — D2é,) 2 ua ,t) = E,G(t) — (4 + és) . (19) 


Then, from Eqs. (2), (15) and (17) 


*) , ; ) 
Kk | Per R = P(r)u an |- = «Rh, - 2 | Po ou — KU £ | Po 2 a> = | 
or or or or or or ar (20) 


So, from Eqs. (5), (6), (15), (18), (19) and (20), one finds 


: fa : a 
(B, — D,é,)(A2 - c( + cd!) | S(r)uR, dr 
ot Jk; 
Cul Le , t) OR Da) 
= «(B, — D,é(As — ¢ ePUL| R ppt? . | Rae) | 
or or 
A ae ED p . ORAL) 
- fe, = DANA, = cat) PUL)| (p,) Sn ® _ en, oRa(E | ' 21) 
‘ or or 
= — K(B, = D f P(L.) OR (L2) @ oo C.&) + (B, = D.£:) 2 Nutt ’ t) 
or or 
y 92 1 42 .2, 9 
— «K(A, — C2; PUL)Ru(Ls)| (A — C,&) + (B, — D&) 2 Jut, i (22) 
AR (Le d G.(0 | 
= | * ees gE, +- KE, Jw ) = 
B, — D,ti)P(L2) —- | (24 tt) (1) — E,xG.(0) 
+ (A, =< 2s \P(L,)R,(L,) Ht “——" Mey Jo(t) aad Ei, xG,(0 . (23) 
Hence 
a } a 
(2 + Ké; | 3, — D,é&)(A2 — Crk) I S(r)uR, dr 
— (B, — D,&,)P(L,) - Oh a (t) — (A, of, P(L,)R.(L,)v 0 | (24) 
€ 
, BRT.) — » " 
= D £, —_ B P( L.) Mattes! EF. xG, t) + (C'£;, —_ A, P\ L,)R,(L,)E, xG; [ 
€ 


This is first order linear, and so may be integrated to give 


wa IRL.) , 
|, - DANA, - CH | S()uR, dr — (B, — D,#)P(L, ) SE) w(t) 
JL, or 


— (A, — C,§,)P(L,)RAL,)vl Is exp (xé,t) 


= (B, — D,&,)(A, — Ck) | S(r)u(r, O)R, 
VL 
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Rs (Ls) 


—3, ~—w0) — (As — ofi)P(L,)R.(L,)v(0) (25) 


— (B, — D,ti)P(L.) 


AR (L2) 4 


+ (D,ki — B,)P(Ls) “EY Bye [ Gi(1) exp (xé21) dr 
Cc J0 


+ (Ce? — A,)P(L,) RAL) Ey « / G,(1) exp (xé2z) dr. 
0 


If this equation is divided through by exp (ké#), the left side will be the same as the right 
side of Eq. (16). Hence 


‘ Ls 4 
T co B, — D,t)(Az — Ck) [ S()Ri dr — P(L.)¥,(B, — D,&) oh, (I) 
J Ly or 
als 
P(L,)b,(Ao — Cok Rate) | = |, — D,ti)(Az — Ck) | S(r)u(r, 0)R, dr 
JLy 
. eS OR,(L.2) y ¢2 ‘ 2 
- (B. — Dt) P(L.) — oe w(0) — (A, — Co€)P(L,)RAL,)v0(0) | exp (— x&d) (26) 
+ (D.; — B,)P(L,) = AE) Ex | G.(r) exp [x&(7 — 0] dr 
or J0 
Cg. — A,)P(L,)RAL,)E,« | G,(r) exp [xti(r — 0] dr. 


This equation may be used to obtain 7',(¢) when R,(r) is known. Its use will be illustrated 
by solving the problem stated in the introduction. 

A formal solution to the problem of heat flow in a cylinder. With the results de- 
veloped above, it will now be relatively easy to solve the heat flow problem stated in the 
introduction, and, incidently, a number of related problems. 

Let u(r, t) be the temperature in the cylinder at a distance r from its axis, and w(t) 
the temperature of the surrounding tube. Then the problem may be formulated math- 


ematically as follows: 


" j (- aw) _ ou . (27) 
rar or dl 

u(r, 0) = f(r), (28) 
w(Q) = W, . (29) 

au(O, 2) 
uO, ) _ (30) 

or 
an ow a[w(t) — u(L, 0], iad 
i 

h dw = alu(L, t) — w(t)] + b[g() — wld). (82) 


dt 
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In these equations k, a, h and g are positive constants, f(r) is the initial temperature of 
the cylinder, W is the initial temperature of the tube, and g(¢) is the temperature of the 
surroundings. Equation (32) arises from physical considerations. However, for math- 
ematical purposes it is preferable to replace it by the following equation, which may be 


obtained algebraically from Eqs. (31) and (32): 


a(b — hé2)u(L, t) + (a + b — het) ee = abg(t) — ai( 4 + cet (0). (33) 

In order to solve for u it is merely necessary to find the radial eigenfunctions, R,, , 
and the equation satisfied by the eigenvalues, &, . The solution can then be completed 
by means of the results developed above. The F, and &, will be obtained by the method 
of separation of variables. To apply this method, one must first set g(t) = 0. The R, 
and £, obtained in this way will still apply for g(t) ¥ 0. 

Let u and w of Eqs. (27)—-(33) be represented by the series given in Eq. (1), where 
u, and w, are assumed to be of the form given in Eq. (15). Then, for g(t) = 0, the u, 
and w, will each satisfy Eqs. (27)—(33), and so by the method of separation of variables, 
one finds from Eqs. (27) and (30) that 


u, = E, exp (—x&)Jo(é.r) when g(t) = 0. (34) 


From Eq. (31), then 


aw, = E, exp (— xé&,d[aJ.(é,.L) — &,.J0,(&,L)] when g(t) = 0. (35) 


Hence, from Eq. (33) with g(t) set equal to zero, one obtains 
a(b — h xé.) J o(é,L) =E(at+ b — heé,)J,(E,L). (36) 


This is the equation which determines the eigenvalues, £, 
Now that the eigenvalues have been determined, the standard separation of variable 
procedure will be discontinued. The restriction that g(t) = 0 will no longer be imposed, 


and Eas. (34) and (35) will be generalized to 
“= Tp SEs), (37) 
aw, = T,(b[aJ(€,L) — &,J,(&,L)]. (38) 


Equations (27), (30) and (33) are special cases of Eqs. (17), (18) and (19) with 


Po) =r, Qn) =0, SH=r, A= C= dD, =E£, = L, = vd = 0, (39) 
A, = b/(hx), B, = (a + b)/(ahx), CG, = i, D, = l/a (40) 
E, = b/(hn), G.(t) = g(b), Le = L. (41) 


With the particular choices for constants and functions given by Eqs. (39)-(41), it will 
now be seen that the u, and w, of Eqs. (37) and (38), with & given by Eq. (36), satisfy 
Eqs. (1)-(6). Moreover, Eq. (31) agrees with Eqs. (4) and (1), and Eqs. (27), (30) and 
(33) with Eqs. (17), (18) and (19). Hence all the conditions used in deriving Eq. (26) 
are now satisfied, and so 7’,(t) may be obtained from Eq. (26). Upon using Eqs. (15), 
(28), (29), (37) and (38), and changing subscripts from k to n, one finds 
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TAO) 
exp (— xé&,t) (42) 


a(b—ht) | rf(r) J (é.r)dr+ah xé,W LJ (EL) +ab xt, LJ (EL) [ g(r)exp( ké&7r)dr 








a(b—hxé?) [ “rJ2Er)dr-th xL [ad o(EqL) — End (EnL) eud sEad) 


Temperatures u and w are now found by combining Eqs. (1), (37), (38) and (42). 
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Models of man—social and rational; mathematical essays on rational human behavior in a 
social setting. By Herbert A. Simon. John Wiley & Sons, Inc., New York, and Chap- 
man & Hall, Ltd., London, 1957. xiv + 287 pp. $5.00. 

The book consists of a collection of sixteen articles, previously published in various statistics and 
social science journals, together with introductory material which has been added to clarify the articles’ 
mutual relations and unity. They are grouped under four main headings: I. Causation and Influence 
Relations; II. Social Processes; III. Motivation: Inducements and Contributions; IV. Rationality and 
Administrative Decision Making. The collection represents a concentrated effort to bring to bear on 
various social science theories the powerful symbolism of mathematics, especially mathematics of a 
non-statistical nature, e.g. algebra and calculus. In view of the considerable range of both subject matter 
studied and methods employed, the book should be very interesting to anyone concerned with the in- 
vestigation of social science topics by means of mathematics. 
G. MorGan 


Light scattering by small particles. By H. C. Van De Hulst. John Wiley & Sons, Inc., 
New York, and Chapman & Hall, Ltd., London, 1957. xiii + 470 pp. $12.00. 


This book is a treatise on the scattering of light in which relatively few phenomena are dealt 
with. The assumptions made throughout are that the frequency is unaltered by the scattering process, 
that ‘the scatterers are independent and without any cooperative effect among them, that there is no 
multiple scattering present. This book is concerned mainly with the scattering theory for a single par- 
ticle and it gives first of all general theorems for particles with arbitrary size and shape. It is shown 
that the scattering by any finite particle can be expressed by four amplitude functions which are com- 
plex functions of the incidence and scattering directions, and knowing these functions permits calcu- 
lation of all quantities involved. 

The book contains the amplitude functions and cross sections for various particles, and finally it 


Ls 
deals with specific applications in a number of branches of physics including various problems in chem- 


istry, meteorology, and astronomy. 
Roun TRUELL 


Physical properties of crystals. By J. F. Nye. Oxford University Press, New York, 1957. 
prot : s : 


xv + 322 pp. $8.00. 
The purpose of this book (as stated by the author), is to formulate the physical properties of crys- 
tals systematically in tensor form. The book is an excellent detailed survey of macro crystal physics; 


there is no attempt to deal with the quantum theory of solids. 

The chapter headings are as follows: (1) The Groundwork of Crystal Physics, (2) Transformations 
and Second Rank Tensors, (3) Para and Diamagnetic Susceptibility, (4) Electric Polarization, (5) The 
Stress Tensor, (6) The Strain Tensor and Thermal Expansion, (7) Piezoelectricity. Third Rank Tensors, 
(8) Elasticity. Fourth Rank Tensors, (9) The Matrix Method, (10) Thermodynamics of Equilibrium 
Properties of Crystals, (11) Thermal and Electrical Conductivity, (12) Thermoelectricity, (13) Natural 
and Artificial Double Refraction. Second Order Effects, (14) Optical Activity, plus seven appendices 
on various topics related to the text. 

This textbook is recommended to both physicists and mathematicians interested in crystal physics 


problems. 
Roun TRUELL 
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SOME STATISTICAL PROBLEMS ENCOUNTERED IN A THEORY 
OF PINNING AND BREAK AWAY OF DISLOCATIONS* 


BY 
G. F. NEWELL 


Brown University 


1. Introduction. From a physical model proposed by Granato and Liicke [1] to 
describe the mechanical hysteresis in crystals caused by the motion of pinned dislocation 
lines, the following type of statistical problems arise. One distributes points (impurity 
atoms) at random along a line (edge dislocation) according to a Poisson distribution 
and then extracts from this line, segments of length Z, a random variable having a 
probability density P(Z). We wish to calculate: 

(1) the probability that in this segment, no pair of second nearest neighbor points 
are separated by a distance greater than some preassigned distance D; the ends of the 
segments are treated as additional points on the same basis as the randomly distributed 
ones. 

(2) the probability distribution for the distance between consecutive points (again 
including the ends of the segment as additional points) for those configurations of points 
satisfying the condition given in (1). 

For a detailed discussion of the origin and physical significance of this problem in 
relation to the theory of dislocation lines, the reader is referred to the work of Granato 
and Liicke which deals mainly with the physical aspects of hysteresis in solids. Although 
the connection between the above more abstract formulation of the problem and the 
physical model described by them will be outlined briefly, the main emphasis here will 
be on the mathematical solution. We shall obtain solutions in the form of single contour 
integrals involving the Laplace transform of P(Z). In some cases, particularly if P(Z) 
is exponential (Poisson distribution of Z) the solution is obtained exactly in closed form. 

A dislocation line is an imperfection in the crystal. The physical details of its structure 
are not important to the present discussion; suffice it to say that it behaves in most 
respects like an elastic string under tension. Displacement of the dislocation line is 
inhibited by impurity atoms and by intersections with other dislocation lines so that 
under stress the dislocation line is pinned at various points along its length and is de- 
formed much like an elastic string fastened at many points. A distinction is to be made 
between “‘weak pins’’ represented by the impurity atoms and “strong pins’’ represented 
by the crossing dislocation lines. The weak pins are capable of holding the dislocation 
line only if the force exerted on them is less than some maximum value. If the force 
exceeds this maximum, the string breaks away from the pin. The strong pins, on the 
other hand, fasten the dislocation line under arbitrary stress. 

If a dislocation line is placed in a uniform stress field which exerts a certain lateral 
force per unit length of line, this force is transmitted to the pinning points. Each pin 
must support half of the force exerted on the two pieces of dislocation line on either 
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side of it. The force on the pin is therefore proportional to the distance between its two 
nearest neighbor pinning points, the total length of line that this pin must help support. 

The dislocation line breaks away from the weak pin if this force exceeds a certain 
critical value or if the distance between its two neighbors on either side is greater than 
some critical distance D which is inversely proportional to the force per unit length of 
dislocation line. If for any given stress a weak pin should break, its load must be supported 
by the adjacent pins. If the adjacent pins are also weak pins, they also break for they 
must now attempt to support a still longer length of line than the pin that broke. The 
result is that all weak pins on either side will break until the line is finally supported 


by the two nearest strong pins on either side. See Fig. 1. 


(a) 
C-e—e"  ~eeo-e ee 0 


(b) 
Ow eee —e—O 
eee, 


Fic. 1. (a) Shows a dislocation line under small stress with no break away. The circles denote strong 


he dots weak pins and the line is supported at all pins. 


pins, t 
and unstable stage in the break away. The weak pin supporting the largest 


(b) Shows an intermediate 
total le ngth of line was broken. 
(c) Shows the final stage of break away in which the line is supported only by the strong pins for that 


segment in which break away has occurred. 


If, therefore, in a segment between two strong pins, there is any weak pin supporting 
pieces of the dislocation line having total length greater than D, all the weak pins in that 
segment break and the segment is entirely supported by the strong pins. Equivalently 
all weak pins between two strong pins break if any pair of second nearest neighbor 
pins are separated by a distance greater than D. 

If weak pins and strong pins were spaced at regular intervals along the line, then 
for any given stress, we would find that either all weak pins in all segments broke or 
none broke. This would give rise to a sudden change in the physical behavior of the 
system at some critical stress, a situation which is not entirely satisfactory from the 
physical point of view. If, on the other hand, the pinning points were distributed accord- 
ing to some statistical law, there would be a certain non-zero probability that some of 
the weak pins would break for any value of the stress. This type of model leads to a 
hysteresis effect which, though strongly dependent upon the amplitude of the stress 
field, does not vary discontinuously. 

Granato and Liicke suggest that a reasonable assumption would be that weak pins 
are distributed at random along the dislocation line, i.e. according to a Poisson dis- 
tribution. Because of interactions between strong pins, however, they proposed to con- 
sider the distances between strong pins as independent random variables, not necessarily 
of Poisson type. They suggest, in fact, that it might be desirable to consider a network 
of strong pins spaced at equal intervals. We will consider an arbitrary probability 
density P(Z) for the distance Z between strong pins but will find that a Poisson dis- 
tribution for Z leads to the simplest results mathematically. 

The questions that now arise are the following: 

1.What is the probability that in any segment of line between two strong pins, no 





to 
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break away has occurred for some given stress or equivalently for some given value of 
the critical distance D? 

2. What is the probability distribution for the distance between consecutive un- 
broken pins for any given value of D? 
These two problems are the same problems as those described at the beginning of this 
section except they are stated in a different form. The second problem above can be 
decomposed into two parts; to find the distribution of distances if break away has not 
occurred and that if break away has occurred. The former is equivalent to the previous 
problem 2 whereas the answer to the latter follows directly from the fact that the distance 
in question is Z having a given probability density P(Z). 

For the subsequent calculations, we introduce the following notation (see Fig. 2): 





Oo 1 2 n-l on ne! 

ro | 42 | | yn | Yrs 

| Z *| 
Fic. 2. In a segment containing n weak pins, points are numbered consecutively from 0 to n + 1 with 
the ends (strong pins) numbered 0 and n + 1. Distances between points are denoted by z;,..., 2n 41 
and the total length is given by Z. 


is a random variable designating the number of points (weak pins) within some 


Ht 
segment of length Z (between strong pins); 
denotes the distance between consecutive points with z, the distance from the end 
of the segment to the first point and z,,, the distance from the nth point to the 
other end of the segment; 
| is the average distance between randomly distributed points; /~* is the density of 
weak pl ~ 
From the definition of z; , it follows that 
n+1 
Z = z Z; (1) 
i=l 


and from the assumption that the weak pins are distributed according to a Poisson 
distribution with density I~’, it follows that I-* exp (— 2l~')dz is the a priori probability 
that the distance z between any two consecutive weak pins lies in the interval between 
z and z + dz. The probability, for fixed Z, that there are n points in Z and 2, has a value 
between z, and z, + dz, , 2» has a value between z, and z, + dz, etc. and there are no 


points in the last interval of length z,., = Z — > 2; ,18 

l exp — — ra ) dz, ‘+. dz, = l ; exp (—I7'Z) I] dz, * for 7? > 0 
i=! 

exp (—1°'z,) = exp (—I°'Z), for n = 0. 


Since P(Z) is the preassigned probability density of Z = }>**! z; , we obtain the 


following for the probability that the interval is of length Z with z, between z, and 
z, + dz, , --* » 2+, between z,., and z,., + 2,4; : 


l"P(Z) exp (—ZI") dz, dz, +++ dzn.; , (n > 0). 
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The prescription for evaluating the average value (f) of any function f(n, 2, , 2. , 
2,41) of the distribution of pins is to compute 


f) = , 3 | dz, | dz. ++ | dz,.,P(Z)l" exp (—ZI)f(n, 2; , +++ 5 Zn+1)- (2) 

It will usually be impossible to evaluate multiple integrals of this type unless f is of a 

special form. Fortunately, the quantities of interest here can be represented by products 
of relatively simple factors or sums of such products. 

2. The probability of no break away. The first problem is to find the probability p 

that no pair of second nearest neighbors in some segment are separated by a distance 

1,2,---,n. This is equivalent 


greater than D, i.e. the probability that z; + z;., < D,7 = 1,2, 
D for all j 


to finding the average of a function f which has the value 1 if z; + 2z;,, 
and 0 if z; + z;., > D for any j. Such a function can be written in the form 
g{(@1 + 22)/D)g[(z@2 + 23)/D] +--+ gll@n + 2n+1)/D], (3) 
in which the function g is defined by 
| for 2< i 
g(x) = == (4) 
(O tor o> 1 


In view of the form of (3), it is convenient to measure z in units of D by writing 


zx; = 2,/D, 
so that substitution of (3) into Eq. (2) for f gives 
p= > D(D/I" | dex, ++ | dit,..P( D > 2,] exp (-—pr > 2) 
n ‘ J0 \ j=1 \ 1=1 


x Il g(a; + 
j=1 


In these integrals, only the function P fails to be in a form that can be factored into a 
product of functions of only one or two integration variables. We can obtain a factored 
form of the integrand at the expense of an extra integration by using the Laplace trans- 


form of P. If we define 


A() = | dZe*P(Z), (6a) 
then 
7 l oe sZ bs 
P(Z) = ar ds e’” A(s), (6b) 
Qxi Ja-io 


in which y is some real number for which A(y) is defined by Eq. (6a). Equation (6b) 
is the well-known formula for inverting the Laplace transform. 
If we substitute Eq. (6b) for P(Z) in Eq. (5), let 


a =(s—l1)D, (7) 


and interchange some integrations and summations, we obtain 





to 
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p = 1 [ da AaD"' +1) >> (/y" [ dz, e*** 
: 0 


n=0 


(8) 


i) 


x | dx,e°**g(a, + x) ++ [ dx, e**"g(2n-1 + Zn) [ d254, 0°" gz, + Sas); 
“0 0 0 


in which y must now be chosen so that A(I~-' + yD~’) exists as defined by Eq. (6a). 

The multiple integrals of Eq. (8) are in such a form that most of the integrals represent 
only a repetition of the single integral operation. Equation (8) can be simplified, formally 
at least, by introducing an integral operator K defined by the equation 


Kf(x) = | dx’ e*' g(x + x')f(x’), (9) 
for any function f(x). In terms of this operator K, Eq. (8) takes the form 
I = / —1\ ” az = n n . 
p= 55 [. daA(aD™ + 1 » | dx e** D (D/IRT. (10) 
K"I denotes » repetitions of the operator K on the function f(z) = 1. 


The sum over 7 in Eq. (10) can be formally written as 
> (D/)"K" = 1 — KD/I". (11) 


Thus the problem of evaluating the sum of multiple integrals in Eq. (8), reduces to an 
investigation of the properties of the operator (I — KD/I)™. 

The operator (I — KD/1) is well defined by Eq. (9) and transforms any function 
v(x) upon which it operates into a function g(x) according to the rule, 


¢ I — KD/)Y(2) = x) — (D D | dx’ e** g(x + x')¥(x’). (12) 


The operator (I — KD 1)~’ must describe the inverse transformation 

v(x) = (I — KD/D“'¢(2). (13) 
In other words Eq. (11) applied to any function g(x) produces the solution ¥(x) of the 
integral equation, Eq. (12) We see that for Eq. (10) we wish to know the function (x) 
resulting from g(x) = 1. For future reference, however, we shall solve Eq. (12) for 
arbitrary ¢(z 

Substitution of the explicit form for g(a), Eq. (4), into Eq. (12) gives 
g(x) = V(x) for x>1 (14a) 


ni-s 


g(x) = (x) — (D/D) dx’ e** W(x’) for x <1. (14b) 
Equation (i4a\ needs no further study. Equation (14b) can be transformed into a 
familiar type of differential equation by first differentiating with respect to z to obtain, 


WD) WE) 5 Doty — 99. (15a) 
dx dx l 
Replacement of « by 1 — 2 gives 
dgp(l1 — 2) _ _dWi-a, Da. = 








160 G. F. NEWELL [Vol. XVI, No. 2 


By solving Eq. (15a) for Y(1 — 2) and eliminating it from Eq. (15b), one obtains the 
equation 
P(x dy(x) , (2) - d°o(x) dg(x) D , dg(1 — x 
—- +a— “Y(2) = > — eT aoe (16 
dx . ew dx” o dx ih dx » 


This equation to find y in terms of ¢ is the famous linear inhomogeneous differential 
equation with constant coefficients, the formal solution of which is elementary. For 
the boundary conditions, we observe from Eqs. (14b) and (15a) that 

W(1) = ¢(1) 
dy(x) D deg(a (17) 


= ——e*y(1) | 
dx aan l al +r ax 


Actually it is not necessary to find the complete solution of these equations. All 


we need in Eq. (10) is the value of the integral 


| dx ¢ ** )(2) = | dx ¢ V(x) os | dx ¢ V(x) 


v1 


t= [ dx e“*g(a) + [y(O = (0) | D) 


“1 
by virtue of Eq. (14a) and (14b). Thus we need find only ¥(0) evaluated for g(x) = 1. 
The explicit evaluation of the indicated steps is straightforward and leads to the 
result 
] ii ; e° —e + (b — a)(l/D) l e~ 
Dp = da A(aD a - — , (18) 
‘ 21 in a + t be” — ae D a } 


where a and b are the two roots of the equation 


af - -(3) = (3) - (2) | : (19) 


We previously imposed a restriction on y which is equivalent to requiring that the 
path of integration from y — i ~ toy +7 © in Eq. (18) lies to the right of all singulari- 
ties of A(aD~* + I”*). The expression in the bracket of Eq. (18) also has singularities, 
infinitely many of them, and they arise originally from the fact that the sum over n 
in Eq. (8) will not necessarily converge if the real part of a is too large. To justify the 
various formal steps leading to Eq. (18) we must be sure that the series converges for 
all a along the path of integration. This is achieved by picking the path to the left of 
all singularities of the expression in the bracket. It is always possible to pick such a path. 
One can show for example that y = — D/l is always satisfactory. 

At this point, we could express p directly in terms of the given probability distri- 
bution P(Z) by using the definition of A, Eq. (6a), to eliminate A(s) from Eq. (18) in 
favor of P(Z). This would give a double integral over both a and Z which would be 
very difficult to reduce to simpler form. The integral over Z cannot be done unless 
P(Z) is explicitly specified and the integral over a is too complicated to perform first. 
We cannot conveniently proceed further unless we specify the form of P(Z) or equiv- 
alently its transform A(s). 

Despite the complexity of Eq. (18), there is a sizeable class of functions P(Z) or 
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transforms A(s) for which this final integration can be performed explicitly. The 
trick is to consider functions A(s) having only a finite number of poles, to close the 
contour in Eq. (18) with a semi-circle in the left half plane and use residue theory. One 
then avoids the complicated singularities of the integrand in the right half plane. 

The class of functions P(Z) which can be treated in this way include exponentials, 
polynomials in Z times exponentials or any linear combination of such. One could at 
least approximate any P(Z) with functions in this class but the analysis would be ex- 
tremely tedious for any but the simplest examples. 

If we have a Poisson distribution of strong pins, i.e. 


P(Z) = L-*e*, (20) 
with L the average distance between strong pins; L~* the density, then 
A(aD"' + U') = DL ‘la + DI + BT". 


In Eq. (18) there is only one singularity to the left of the contour, namely the simple 


pole of A ata = — DIl-' — DL™. The value of p is found explicitly as 
D fe" —e" +(b— aD) _ wee ‘s 
P L be" — ac F D a ‘ (21a) 
with 
a\ a a\" DY “ si 
-* -(2) + (3) — (2). | ; (21b) 
nd 
a= —D(l'+L"'). (21c) 


The probability p is a function of only two independent variables D/l, the ratio 
of the break away length to the average distance between weak pins and D/L, the ratio 
of the break away length to the average distance between strong pins. The values of 
these two ratios determine a, a, and b as given by Eqs. (21b) and (21c) from which p 
is evaluated. 

\ discussion of these results and some approximations derived from them will be 
described in Sec. 4. 

3. Average of additive functions. Most practical problems are concerned with 
finding the average value of some function which is the sum of contributions from each 
of the intervals of dislocation line between consecutive pins. We therefore desire to 
compute the average (G@) of a function of the form 


G(z,;) + G(ze) + +++ + G(Qn41) (22a) 


| 


G(n, 2; , 22, °** 52) = if z;+2;.,< D for allj 
Laz) if 2z; +21, > D for anyj. (22b) 
This class of functions includes as special cases the following: 
1. If G(z) = z, then both Eq. (22a) and Eq. (22b) have the value Z and (G@) = (Z) = L 
independent of D. 
2. If G(z) = 1, Eqs. (22a) and (22b) give the values n + 1 and 1 respectively. The 


value of (G@) — 1 is therefore the average number of unbroken weak pins per segment. 
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3. If G(z) = 6(z — 2) then (G)/Jf> (@) dz is the probability density for the distance 
between consecutive unbroken pins. The evaluation of this was listed in Sec. 1 as one of 
the problems of particular concern. Knowing this for arbitrary z, we could in fact go 
back and evaluate (G@) for arbitrary G by superposition. We shall treat the problem 
of finding (G) for arbitrary functions directly, however, since the determination of the 
probability density for the distance between unbroken pins is usually of no direct 
physical importance itself but only a means of finding the average of certain additive 
functions. 

Of particular interest to the theory of mechanical hysteresis is the case G(z) = 2’. 
The displacements of dislocation lines contribute to the strain and consequently also 
to the elastic coefficients. For the loading portion of the stress-strain curve, this contri- 
bution to the elastic coefficient is proportional to the sum of cubes of the interval lengths, 
i.e. to (G) with G(z) = 2°. The stress-strain curve for loading is non-linear because (G) 
depends upon the break away distance D which in turn is inversely proportional to the 
stress. As more pins break with increasing strain, the ratio of strain to stress increases. 
The unloading portion of the stress-strain curve is linear, however, because those pins 
that are broken on the loading portion of the curve remain broken during the unloading. 

With the aid of Eq. (4) we can write Eq. (22) in the form 


n+l ’ 
Gn, 2 5 *** 5 Sar) = GZ) — | a — > ce) | I] g{(Ze + 2x+:)/ DI. 


(G@) ean then be split into three parts by writing 


(G) = (G;) + (G. _ G. (23) 
with 
G)= DL] da-++ | dens:P(Z)I™ exp (—Z/D 
n=O JO 79 (23a) 
n+1 r, 
x | = ae) | IT gle. + 2x-:)/D) 
7=1 =1 
G:) = | dZP(Z)G(Z) (23b) 
G)= DO] da | dzn.:P(Z)I exp (—Z/DG(Z) TT gl(ee + zus:)/D]. (28) 
(G,) and (G,) — (G;) give respectively the contributions from pins that do not break and 


those which do break. 

In view of the fact that neither P(Z) or G(Z) is explicitly specified, Eq. (23b) cannot 
be simplified any further. (G;) is easily handled because Eq. (23c) and Eq. (5) differ 
only in that Eq. (23c) has P(Z)G(Z) where Eq. (5) has P(Z). Since P(Z) is arbitrary, 
we can evaluate (G;) simply by replacing P(Z) in the formula for p by P(Z)G(Z) or 
equivalently by replacing the transform A(s), Eq. (6), by 


Bis) = | dZe~"* P(Z)G(Z). (24) 
0 


We can also write (G,) = B(0). 
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Only (G,) remains to be evaluated and this we do by a method which parallels that 
of the last section. In Eq. (23a), the summation over 7 is taken outside the integrals; 
z; is replaced by Dx; ; P(Z) is expressed in terms of A(s); and the integrals are written 
in terms of the operator K noting that for any j, G(z;) appears only in the integral over 
z, , all other integrals being of the same type as in the last section. This leads to the 


Zi ys 


form 


l e 7 - 
G.) = 5— | da A(aD™' + I"') 
Gat Jy-;5 -_ 
(25) 


2 n+1 


x [ dx, nd ye (D/o"K'a(Da 1}. 


n=0 j=1 
The expression K’~'G(Dx)K"~**'I represents transformation of the function I with 
K,n — j + 1 times, multiplication of this function by G(Dz) and the transformation of 


the product with K, 7 — 1 times. 
If we rearrange the order of summations by letting k = 7 — 1 andm =n —j +1, 
the expression in the bracket can be written as 


> ¥ (D/)'K'G(D2)(D/)"K"l = (I — KD/D“'G(Da)(1 — KD/D)". 


The function ¥(x) = (I — KD/l)"'I was encountered in the last section. It was also 
shown there how (I — KD/1)~'g(x) could be found for any function ¢(x). This must be 
done for g(x) = G(Dzx)y(x) as required above and the integral over x, in Eq. (25) evalu- 


ated. This procedure leads to the final result that 


. ~1 fe? — e* + (b — a)l/D ‘I 
r( T 3( j_— : da — —— —-——- — - ——  — 
¢ BO a | | aBlaD’ +1 a D 


(be* — ae’) 


ts [ da A(@D™ + I")(be" — ae’)? [ dx G(D2) (26) 


2m . /J0 
x {(b +e°D/)) exp [}(a — b)z] — (a + e’D/D exp [—}(a — b)z]}?. 


The remaining integrations cannot be done explicitly unless the functions P(Z) and 
G(z) or A(s) and B(s) are specified. There is again an important class of functions P(Z) 
and G(Z) for which these last integrations can be performed analytically. They can be 
done if P(Z) and G(Z) are linear combinations of polynomials times exponentials (the 
same class of functions for which p could be evaluated) and, in particular, if G(z) is 
proportional to 2°, k integer and P(Z) a Poisson distribution. The trick is to close the 
contour of the a integral in the left half plane and use residue theory as done previously 
to obtain Eqs. (21) for p. 

\lthough the evaluation of (@) in this latter case is straightforward, the residue of 
the first integral in Eq. (26) requires k differentiations of the quantity in the bracket 
and the evaluation of the integral over z in the last term of (G) involves integration 
of z’ times exponentials, all of which are elementary operations but unfortunately ones 
which lead to very awkward formulas. We shall not give these formulas here. 

4. Limiting cases. Although we have seen that it is possible to obtain results in 
closed form for many physically interesting properties of this model, the formulas are 
awkward to evaluate and it is highly desirable to obtain simple approximations capable 
of giving at least qualitatively correct results with a minimum of effort. 
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For any particular material, the values of /”’ and L~*, the densities of weak and strong 
pins respectively, would be fixed but D~* which is proportional to the stress would be 
variable. The critical range of D in which break away is most likely to occur will however 
depend upon /~* and L~' which themselves vary widely depending upon the material. 
For a very pure material we may find / > L, whereas for a less pure material we might 
have L > 1. In either of these limiting cases, one can obtain simple approximations. 

I. Pure materials. If / > L, there is a high probability that there will be no weak 
pins between two consecutive strong pins and therefore nothing to break. There is a 
small probability that there is one weak pin in the segment between strong pins and 
this pin would break or not break accordingly as z, + z. = Z is greater or less than D. 
We could find p directly from Eq. (5) by evaluating individual terms for n = 0, 1, 
Since there is a very small probability for having more than one weak pin in a segment 
(n > 1), the first two terms of Eq. (5) would give a fairly good approximation. 

A slightly better and in most cases simpler approximation results if we note that for 
Z < D, no break occurs for any value of n nor can it oecur for any Z if n = 0. We con- 
clude that 

aD . 


p> | P(Z)dZ+ | P(Z) exp(—Z/) dZ, (27) 


a 


v0 /D 


in which the first term represents the probability that Z < D and the second term the 
probability that if Z > D then n = 0. We also can obtain a close upper bound by noting 
that the probability of break away is at least as great as the probability of break away 
forn = 1, i.e. 


p<l1l- j P(Z) exp (—Z/1)(Z/l) dZ. (28) 


“D 
These integrals can be evaluated for a wide variety of functions P(Z). For example, 


if P(Z) is Poisson then 


; > 1 — exp(—D/L) + (1+ L/l)” exp [-—(D/D)0 + L/D] 27a) 
<1—(L/D0 + L/D [1 + (D/L) + L/D] exp [—(D La + L/)). 28a) 


If P(Z) = 6(Z — L), equal spacing of strong pins, then 


] ] for D> L 27b) 
> ae ee 


f 


1 — (L/I) exp (-—L/D) exp(—L/l) for D< L. 28b) 


These bounds are valid for arbitrary values of L/l but are close bounds (error of order 
L’/l’) only for L/l « 1. For large stress, D — 0, and L/l « 1, Eqs. (27) and (28) show 
that p — 1 ~ L/I for any P(Z). The probability of break away, 1 — p, is therefore 
small, of order L/l, because of the shortage of weak pins, but most of the break away 
that does occur takes place for D of order L. The amount of break away is sensitive 
to 1/L but the stress at which it does occur is nearly independent of / for Ll < 1 


For L/l = 1, Fig. 3 shows the upper and lower bounds of (1 — p) given by Eqs. 
(27a, b) and Eqs. (28a, b), also the exact value of 1 — p for the Poisson distribution of 


Z given by Eq. (21). Even for L/l = 1, Eq. (27a), the upper bound on 1 — p, is a reason- 
able approximation. The accuracy of Eq. (27a) is, of course, better for smaller values 


Ww 
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Fic. 3. The probability of break away 1 — p is plotted as a function of L/D (the stress) for L 1 = 1. 
sroken lines 1 and 2 are the upper and lower bounds on 1 — p for equal spacing of strong pins, Eqs. 
27b) and (28b). Curves 3, 5 and 4 are upper bounds, lower bounds and the exact solution respectively 
for 1 — p with an exponential distribution of strong pins Eqs. (27a), (28a) and (21). 

















of L/l but even for large values of L/l, Eq. (27a) is accurate in the limit of small D 


(large L/D). 

The value of (@) can be approximated using similar arguments. It is clear, first of 
all, that if there are very few weak pins, L/l — 0, or if all the weak pins have broken, 
D — 0, then the distance between unbroken pins is always Z and 


G — | dZ P(Z)G(Z) for L/l-0 or D-O. (29) 


If on the other hand, .D — © (zero stress), nothing breaks and one can find (G) as a 
limiting case of Eq. (23) or Eq. (26). The result is 


“0 


: aZ 
G) - | dZ P(Z) {a2 wae ix | dz G(z)(Z — z + 2be™ 7 (30) 
for Do. 


If, in particular, we let P(Z) = L™' exp (— Z/L), Eq. (29) gives directly 


«EL | dZ G(Z)e"*, L/i-0 or D-O (29a) 


70 


whereas Eq. (30) reduces to 
GG) L"(1 + L/p’* | dZ G(Z) exp |-—(Z/L)\(1 + L/D] for Do @. (30a) 


The similarity of Eqs. (29a) and (30a) results from the fact that if weak pins and strong 
pins separately have Poisson distributions, they combine to give a Poisson distribution 


of all pins with average density L™'(1 + 1/L). 


As a further specialization we find from Eqs. (29a), (30a) that for G(Z) = 2" k 
integer, 


(GG) k!IL" for L/l-0 or D-O, (29b) 
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G) > KIL + L/)'* for Doo. (30b) 


In each case, Eqs. (30), (30a), (30b) reduce to Eqs. (29), (29a), (29b) respectively if 


L/l— 0. 

To determine how (G) depends upon D for L/l < 1 we may start from Eq. (29) and 
make a correction for the presence of a small number of unbroken weak pins or we can 
start from Eq. (30) and make a correction for a small number of broken pins. In either 
case we shall approximate the effects of unbroken or broken pins respectively by con- 
sidering only the contributions for n = 1, one weak pin between strong pins and neglect 


= 


effects of n > 1. This leads to the relations 
»D aZ 
(G) ~ (@)p-n + | dZ P(Z)e" 71" | dx[2G(x) — G(Z)] (31) 
. ae aZ 
(G) ~ (G)p.2 — | dZ P(Z)e"?7"T"" | dx[2G(x) — G(Z)]. (32) 
/“D J0 


Both of these formulae are correct to order L/l for all D but Eq. (31) approaches the 
exact limit given by Eq. (29), for D — 0 whereas Eq. (32) approaches the exact limit 
given by Eq. (30) for D— ~. If P(Z) = 6(Z — L), (G) has a discontinuity at D = L 
as did p. For the Poisson distribution, there is a more gradual change in the range D 
of order L. 

As a special case, we obtain 


(G) ~ kIL'§1 + (1 — kY(L/DA + L/D ik + II! | dx e*x’ f (31a) 


(32a) 


5 


G) ~ kKIL*(1 + L/D'"51 — 11 — A(L/DA + L/D Ok + DI" 


. | dx e~*x"*" > 


JD(i-*+L~*) } 


~@ 


if P(Z) is Poisson and G(Z) = Z”. The integrals above can either be evaluated directly 
or from tables of incomplete I'-functions. 

II. Impure materials. For 1 < L, it is clear that if D is of order L, break away is 
nearly impossible because a negligible number of intervals between weak pins of average 
length 7 will be of length comparable with D as required for a break. On the other hand, 
if D is of order 1, break away is almost certain because there is a high probability that 
at least one of the many second neighbor distances will be greater than D. We conclude 
that the range of D where break away is most likely to occur satisfies the condition 


LKD<«L. (33) 


By neglecting terms of order exp (— D/l) compared with 1 in the equations of Sec. 
2, one finds that for an exponential distribution of Z, Eq. (21) reduces to 


-1 
p~ [ a ; (P 1 )e we (34) 


and for most distributions P(Z) one can use an approximation 


p~ | dZ P(Z) exp  - Z(P = i) on (35) 
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If all segments are of length Z = L, then 


p ~ exp | - ‘ (2 - 1)?" | (36) 


which means that for L/l > 1, p changes rapidly from 0 to 1 when D/1 is approximately 
of order log (L/l). 

One can show from the equations of Sec. 2 or from physical arguments that for 
L/l> 1, break away does not appreciably distort the shape of the distribution of lengths 
between unbroken pins. The reason is that break away results mainly from some ex- 
ceptionally long interval between a pair of weak pins. The existence of one such long 
interval does not significantly distort the distribution of lengths of the remaining intervals 
in Z. Yet when one pin breaks, they all break. The distribution of lengths between pins 
that break is therefore nearly the same as the distribution of lengths between pins that 
do not break. 

To find (@) approximately we simply take (G) for the interval between two unbroken 
pins and multiply it by the average number of such intervals between strong pins if 
there is no break. We add to this (G) for the broken intervals to obtain 


a ia if dz G(z)l el f dZ (1 + Z/l)P(Z) exp E Z(D - 1} 
ins #0 (37) 
x / » 
+ [ dZ PYG), — exp |- Z(D ae 1) ]}. 


These integrals are easily evaluated for a wide variety of possible G(z) and P(z). 

To illustrate the approximations of this section, we have plotted in Fig. 4, 1 — p vs. 
L/D (break away vs. stress) for P(Z) an exponential distribution with L/l = 100. 

5. Conclusions. We have found in Secs. 2 and 3 the exact solution of the problem 
described in Sec. 1 and obtained several limiting approximations. This problein is of 
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Fig. 4. The probability of break away, 1 — p, is plotted as a function of L/D (the stress) for L/l = 
100 and an exponential distribution of strong pins. Curve 1 is the upper bound, Eq. (27a) [the corres- 
ponding lower bound Eq. (28a) is too small to draw on this scale]. Curve 2 is the exact solution, Eq. 
(21), and curve 3 is the approximation for L/l > 1, Eq. (34). Curve 4 (broken line) is the corresponding 
approximation for equal spacing of strong pins, Eq. (36). 
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some practical importance to a field of physics that is still in a rather infant and uncertain 
stage of development but it also suggests a number of interesting mathematical problems 
that are somewhat different from the classic problems of probability theory. 

None of the individual mathematical steps in the above analysis is new although 
they are borrowed from rather diverse origins. The analysis is formally quite similar 
to classic methods for dealing with Markov chains. If P(Z) is exponential, the distri- 
bution of distances between second neighbor points can, in fact, be described in the 
form of a Markov chain. It would have been possible also to carry through most of 
the formal steps of Secs. 2 and 3 even if the z; were distributed according to a Markov 
chain rather than a Poisson distribution. One could, however, not expect to obtain 
explicit solutions in such cases. 

Acknowledgments. The author is indebted to Professor Kurt Liicke and Dr. 
Andrew Granato for bringing this problem to his attention and for numerous very 
enlightening discussions particularly regarding the physical basis of various theories of 


dislocation line behavior and related topics. 
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—NOTES— 


EXPERIMENTS IN THE SMOOTHING OF DATA* 


3y MARK LOTKIN (Anco Research and Advanced Development Division) 


1. Introduction. When using sequences of observed data for the evaluation of certain 
desired functions of these data it is frequently necessary to apply some smoothing 
process in order to avoid unreasonable results. In such cases the points of paramount 
concern are then not only what smoothing process to choose, but also when to apply 
it. The following remarks deal mainly with the “when” rather than the “how.” In 
particular, the question of argument smoothing versus function smoothing is discussed 
here, by means of a few simple examples. Both functional values and derivatives are 
considered. 

While the results obtained may not be of the greatest generality they are, however, 
considered sufficiently indicative to be of some practical utility. 

2. Definitions and specifications. Let the exact values of a function u — f(t) fora 
sequence ¢; , = 1, 2, --- , m, of the independent variable ¢ be denoted by u, , and let 
us assume that there are also available observed values of u, denoted by u’ . Of interest 
are the magnitudes of the errors e,;(u’), e;(u’) = uf — u,; , as well as the length |! e(u’) |! of 
the error vector: 

{|} e(u’) J? = >> ew’). (1) 
i=l 

A process which replaces the u/ by values u* that differ as little as possible from the 
corresponding u’ while diminishing the vector || e(u*) ||, e;(u*) = u* — u; , will be called 
a “‘smoothing”’ process. The particular smoothing process to be considered here is based 
on a moving least squares polynomial of fixed degree [1]: 


p(t) = ao tat+--- +a,t. 
The smoothing consists in putting u* = p(t,;), U* = (dp/dt), . For the sake of simplicity 
we assume that m = 2n + 1 is odd, and that the sequence ¢; is equidistant: t; = ¢, + 
(¢ — 1)At,2 = 1,2, --- , 2n +1. 
Then the usual transformations 
i nae (t; ioe tr+i)/At, 
Vi-n-1 = Uj, $=1,2,---,n+l, 


move the midpoint from ¢,,; to tas; = 0, change the ¢; tor, = ¢ — n — 1, and center 
the uf about », . Further, p(t) transforms into a polynomial r(r), so that 


p(t.) = (7) _ PB QT) ’ 
p=0 
and, consequently, 
ux, = 2(0) =a, (2) 
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U*., = a,/At. (3) 


The coefficients a, , p = 0, 1, --- , r of the least squares polynomial x(r) are obtained 
in the well known manner from the system of conditional equations Ca = v, by solving 


the associated normal equations Na = w, with N = C’C, w = C’v, where the super- 
script 7 denotes transposition. Thus, if we put 
Ao = av, a, = adv, (4) 
it is seen that, for example, for r = 3, n = 2, 
a, = (1/35)(—3, 12, 17, 12, —3), 


a, = (1/12)(+1, —8, 0, +8, —1). 


Let us now consider the following special smoothing problem. Suppose there are 
‘ of observed values, and a function f; = f(z/, y’). In applying 


given two sequences 2} , 7 
a number of possibilities present them- 


smoothing procedures to the caleulation of f 
selves: 
I. Function smoothing. Smooth the functional values f! only. 
II. Argument smoothimg. Smooth the arguments 2% , y/ , into «* , y* then calculate 


f*¥ = f(x? , y%). 


/ 


III. Argument and function smoothing. Smooth the x/ , y/ , and subsequently calculate 


and smooth also the f* . 
Since the amount of computational labor involved in each case may vary greatly, 
appropriate indications as to which of these possibilities may be the most economical 
would be of considerable value. While the experiments described below do not provide 
a general answer to this question, they do permit certain interesting conclusions. 

3. Outline of procedure. In order to investigate this problem we conducted the 
following experiments. First exact mathematical functions z(t), y(t), f(0) f(x(t), y(d) 


were chosen, and their values, as well as those of X(¢ dx/dt, Y(t), F(t) were computed 


for a sequence ¢, of equispaced values of ¢. 
Next the x; and y; were subjected to certain well defined, preferably pseudo-random 


perturbations which produced corresponding sequences zx’, y’, f’ = f(x’, y’). 

n using Method I, Eq. (4) was applied tov = f’, so that the a produced smoothed 
quantities (f’)* g of the function. On the other hand, in using Method II, Eq. (4a) 
was applied in succession tov = 2’, v = y’, resulting in smoothed arguments 2*, y* 
From these there was subsequently calculated f* = f(2*, y*). In applying ITI, finally, 
the smoothing process was imposed also on v = f*, now leading to “doubly”? smoothed 
values (f*)* = h of the function. 


In computing derivatives by means of a, of Eq. (4) the smoothing of f’ led to the set 
g of derivatives of f’, while the smoothing of f* resulted in quantities y of derivatives of 
f*. A comparison of items ¢, ¥ with F then indicated the relative accuracies obtained. 

The computation of rates of change of f may obviously be accomplished also in a 
number of different ways. Thus, for example, quantities X¥*, Y* may be obtained first, 
by the application of the process a, to x*, y*, and then a derivative of f calculated by 


means ot 


: ; 9 ; 
Ft = (z* y*)X* + Of (* ys) Y*, 
: : oY an 
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4. Examples. A. Example 1. We took 
x(t) = sin (t/10), y(t) = (1/10) exp (4/10) 
fz,y=at+y, 
fort = 0,1, --+ , 30. 


Now x and y were perturbed as follows. Denote the third digit of x by d, . Then 
we put 


t’ = (t/10) + e(—1)""d; , «= .001, 
2’ = asin i’, y’ = (1/10) exp ?¢’. 
ff=2'’+y’. 
The calculations showed that, for 2 < t < 28, 


e(t’) ||? = 915.10°° 


e(x’) ||? = 450.10°°, | e(y’) ||? = 520.10°° 
e(f’) ||? = 480.107°. 
Smoothing of t’ by Eq. (4) resulted in |! e(¢*) |!" = 395.10°°, so that the smoothing effect 
on t’ itself was very pronounced indeed. 
Smoothing of x’, y’, f’ into sequences 2*, y*, f* = 2* + »* had the following effect: 
e(x*) ||? = 170.10°°, || e(y*) ||? = 130.10°°, | e( f*) ||? = 240.10°°. 


l 


Since f = x + y isa linear function, g = f*. Further smoothing of f* led to h, and it was 
found that || e(h) ||? = 160.10°°, to be compared with a value of |! e(f*) ||? = 229.10°° 
for the corresponding entries. For the derivatives we found 


e(g) ||? = 326.10°°, | e(y) ||? = 154.107°. 


Due to the linearity of f, F* = y. 
B. Example 2. Next, we let 


a(t) = (1/2)(t/10)’, y(t) = (1/10)(t/10)* + (t/10) 
(Oo =y- 72, 


fort = 0,1, --- , 30. For the disturbances £, 7 of x, y, respectively, we chose the following 
functions: 


ae=a,t+8,, Y=HY¥+n, 


€(-1)"'o, a = (-D™ 1, , 


Il 


wtr 


with « = .001, and a;, , r;, denoting, respectively, the seventh decimal of sin [(60 + 2) 
x/180], tan [(50 + 72) 2/180], forz = 0, 1, 2, --- , 30. 
It was seen that 


e(x’) ||? = 715.10°°, | e(y’) ||? = 708.10°°, e(f’) ||? = 1862.10~° 


between 7 = 2 andi = 28. 
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Smoothing by Eq. (4), 


For 4 <2 < 26, e(j* 
computations led to 


e(w) 


*) [[* = 358.10, 


NOTES 


leading woz. 


|| e(y*) 


~ = 938.10", 


1? = 437.10°° 
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and f* had this effect: 


= 468.10°°, || e( f*) ||? = 1428.10°°. 


C. Example 3. Next we considered the following example: 


a(t) = (t-— 1 t+ ] 


to be 
manner: Let o4, 75 


with « 


andi = 28, 


e(x’) 11° = 375.10 


Smoothing of x’ and 7’ led to 2%, 


e(a*) 


) = At+1 


Y\C) ~ 


( 
2 = a= 


| > 
om 


0005. Similarly for 1 It turned out that for 2’, 1 


, 2 


Evy | 


= 284.10™," 


| e(h) ||> = 864.10°°. As to derivatives, the 
|| e(g) ||? = 1351.107°. 
(2 + 1), f(x, y) 
zey = 2t — 1)/2t + 1), 


omputed for t = 3, 4, --- , 30. The “roughing up” of x, y was done in the following 
denote the fourth decimal of z, y, respectively. Then we put 


for o, #0 


for 4a; 0, 


x’y’, between 7 = 5 


525.10° 


y*, and f*. It was seen that 


I| e(y*) ||? = 420.10°°, 


so that the smoothing did diminish somewhat the length of the error vectors. Further, 


e(f*) ||? = 526.10"° 


The values of g, obtained by smoothing the functional values f’ showed very close 
agreement with those of f*. Further smoothing of f* resulted in values h, for which 


e(h) |\° = 268.10 


*. which was to be compared with 305.10~° for the corresponding 


entries in f*. The derivatives of f, computed in the three different ways described in Sec. 


3. led to 


e(W) 236.10°° 


D. Example 4. Finally, 


r( ft 

y(t) = 

f(t 
fort = — 15, — 14,°---,+ 
“<inusoidal’’ manner: 


< — 
.( 


forz =~ 0,+1,+2,::: 


- (1/10)((t 


e sin (7i/3)/sin (r/3), 


; and E 


} e\¢g) 


we took 


10 


—(1/10)((¢ 


= a(t)/y(t), 


15. A perturbation was introduced into z, y in the following 


UE 


— 3)7((t 


10 


540.10 | e(F*) ||? = 248.10~ 


10) + 1)° 


— 3)°*((t/10) + 1) : 


Yi = Yi +n, 


= esin [(7/3)(¢ + 1)]/sin (1/3), 


.0005. This resulted in the perturbations 


, 0, —e, —¢, 0, °° 
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fori = 0, 1, 2, --- . We found that, for — 13, < ¢ < 13, 
| e(x’) = 450.107*, {| e(y’) ||? = 450.10, || e(f’) ||? = 386.107°. 
In roughing up x’, y’, f’ at? = — 10, as well as in smoothing these quantities, the values 


had been left unchanged at zero, in order to avoid absurd entries there. 
The smoothing again did diminish the error vectors somewhat: 


e(z*) ||" = 383.10, '! e(y*) ||? = 379.10, || e( #*) ||? = 341.107°. 
On the other hand, e(g) = 370.10. Further, | e(h) |? = 99.10°*, while || e(f*) |? = 
128.10°° for the corresponding entries between 7 = — 11 and7z = 11. For the derivatives 


the following results were obtained: 
| e(y 118.10°“, e(g) ||? = 161.10°°, || e(F*) ||? = 132.107°. 


5. Conclusion. In view of the results obtained above the following conclusions, at 
least for the examples discussed, seem to be justified: 

(a) f* is not significantly better than g; 

(b) the improvement A due to the further smoothing of f* may be worthwhile; 

(c) F* is about as good as y, while ¢ is considerably worse than either of the foregoing. 
In summary, then, it may be stated that if functional values f alone are desired, possibility 
[ of Sec. 2 would seem to be adequate. However, if rates of change df/dt also are of 
interest, then possibility IT should be utilized. 
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TRANSVERSE VELOCITIES IN FULLY DEVELOPED FLOWS* 
By STEPHEN H. MASLEN 
By STEPHEN H. MASLEN (Lewis Flight Propulsion Laboratory, Cleveland, Ohio) 


We consider an incompressible flow in a channel whose generators are parallel and 
lie in the x-direction. The driving force, either a pressure gradient or some kind of a 
body force, is also in the z-direction. If the driving force and the flow velocities are 
independent of .r, then the flow is said to be fully developed. In such a case, it is usually 
assumed that the velocity components normal to the z-direction vanish identically. 
However, in some cases it is not clear that this must be true. One such example is that 
of a free convective flow where the velocity (x-direction) shows a cellular structure [1], 
being sometimes in one direction (up, say) and sometimes in the other (down). It is the 
purpose of the present note to show that regardless of these circumstances the transverse 
velocity components in a fully developed incompressible flow must vanish everywhere. 

Consider a fully developed incompressible flow, where the velocity components 
u, Vv, w, are in the (cartesian) x, y, z directions. If P, p, v, and F are pressure, density, 
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kinematic viscosity, and body force, respectively, the equations become 
, : ’ J» 1 
v,+w, = 0 | 
vu, + wu, + P,/p = vAu+ F 
(1) 
vv, + we, + P,/p = vdv 
vw, + ww. + P,/p = vAw 


where subscripts denote partial differentiation and A is the Laplacian operator. The 


boundary conditions are that u = v = w = O at the walls. 
The first of Eqs. (1) implies the existence of a stream function y, such that v = y, , 
w = — y,. Using this definition, the last two of Eqs. (1) can be combined to read 
L(y) = vAAy + (¥, Ay. — y.Ay,) = 0. (2) 
The appropriate boundary conditions are that y, = y, = 0 at the wall. One can observe 
that 


YL(y) = v(Ay)* = ivvAy, + WyWer = v.V,:) 7% vy,Ay a VAY, T Vi) $s 
+ {Woay, — vit + Ww.) —ry AV tev. (+ WI, 


(wP 
st — wy.Ay — WAVE + ¥ 
p 


\ 


+ 

= 

s. 
+ 

. 


+ {—-=+ — wy, Ay 


If this equation is integrated over the cross section, D, of the channel, there follows, 


using Green’s theorem 


ll [yL(y) — v(Ay)"] dy dz = | | -#2: + vp. Ay + 2¥,(y, + | dy 


+ | -¥ — vp, Ay + 4¥.(¥3 + | dz, 


where C is the boundary. If Eqs. (2) are noted, this becomes simply 


, ir ir 
y | (Ay)* dy dz = +- | W(P, dy + P, dz) = +- VP, ds, 
JD PJc P Jc 
where s is the distance along the boundary in the y — z plane. 


In general the cross section, D, is not simply connected, but has a boundary con- 
sisting of several closed curves c; . Because y, , ¥, vanish on each, it follows that y is a 


constant, y; , on each. Then 


y [ (Ay)? dy dz = ; » y, [ P, ds. 


“D 
Because the pressure must be single valued, the right hand integrals must vanish. 
Hence, so does the left hand one. This is only possible if AY = O everywhere. Then, 
subject to the boundary conditions of Eqs. (2), ¥ can be at most a constant and the 


transverse velocities must vanish everywhere. 
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Hence, in fully developed incompressible flows, there can be no transverse velocity 
components. The flow is then found from the second of Eqs. (1), dropping the inertia 
terms. Of course, for free convective flow (for example, Ref. [2]), an energy equation 
must also be included. 
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ON AN INEQUALITY OF LIAPOUNOFF* 
By AUREL WINTNER (The Johns Hopkins University) 


1. Let p(t) be a continuous function which is positive on the ¢-interval under con- 
sideration. Instead of p(t) > 0, it will be sufficient to assume that p(t) = 0, provided 
that p(t) > 0 holds on a dense t-set. The role of this proviso will be that of excluding 
the existence of a function x(t) which satisfies the differential equation 


x’ + p(x = 0 (1) 


and is a non-vanishing constant on some ¢-interval. If such a solution x(t) of (1) is dis- 
regarded [and (1) cannot have two, linearly independent, such solutions in any case], 
then, besides the continuity of p(t), only the assumption p(t) 2 O will be needed. Only 
real-valued solutions x(t) will be considered, and the trivial solution, z(¢#) = 0, will be 
excluded. 

It is clear from (1) that x(t) = 0 or x(t) S O at a given ¢ according as x(t) < 0 
or x(t) > Oat that ¢. Hence the graph of x = x(t) must turn its concavity toward the 
f-axis at every ¢. Since the clustering of zeros of the derivative 2’(t) has been excluded, 
it follows that the zeros of x(t) separate, and are separated by, the zeros of x’(t) [provided 
that either x(¢) or x(t) has at least two zeros]. Let a closed ¢-interval [c, d] be called a 
primitive interval of x(¢) if neither x(é) nor .x’(t) has any zero in the interior of [e, d]. Such 
an interval [c, d] will be called a complete primitive interval of x(t) if for no «€ > 0 is 
fe — e,d + «| a primitive interval of x(é), that is, if x(t)  Oand 2’(t)  Ofore <t<d 
but either z(c) = 0 and x’(d) = O or x’(c) = 0 and x(d) = 0. Note that z’(c) # 0 and 
a(d) =~ Oin the first case, and that x(c) + 0 and x’(d) # 0 in the second case, since the 
simultaneous vanishing of x(¢) and 2’(¢) leads to the trivial solution. 

2. The purpose of this note is to show that, owing to the concept of a primitive 
interval, a theorem of Liapounoff (see below) can be extended from his ‘‘disconjugate”’ 
case to the general case, as follows: 

If p(t) ts continuous and non-negative on a closed t-interval [a, b], and if [a, b] consists 
of exactly n primitive intervals of some solution x(t) of (1), then 


/ " y(t) dt > n2/(b — a) (2) 
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(here n is a preassigned positive integer, and n = 1 is allowed; but n = © is not allowed, 
since solutions x(t) which are constant on some f-interval are excluded). 

Accordingly, if P is defined by 


s 


P a (0 — a) | p(t) dt? (3) 
and if [P] is the greatest integer not exceeding P, then a solution x = x(t) of (1) on 


[a, b] cannot consist of more than [P] consecutive stretches on each of which both x(t) 
and x’(z) are strictly monotone (nowhere constant). The italicized lemma contains how- 
ever somewhat more, since it does not exclude such linear stretches of the graph x = x(t) 
as are not parallel to the t-axis. The proof proceeds as follows: 

3. First, it is sufficient to prove that (2) must hold when [a, }] consists of exactly 
n complete primitive intervals (rather than of exactly n arbitrary primitive intervals) 
of some solution z(t) of (1). In fact, if the value of n is retained but [a, 6] is replaced 
by [a*, b*], where a* < a < b S b*, then the value of the quotient on the right of (2) 


is decreased, whereas the value of the integral on the left of (2) is not decreased, since 


p(t) = 0. 

Next, it is sufficient to prove (2) for the particular case n = 1. For suppose that 
(2) is true for nm = 1. Then, if [a, b] consists of n primitive intervals, and if the latter are 
denoted by [e;,-,; , ¢.], where c = a,c, = band k = 1, --- , n, then an application of 
the case n = 1 of (2) to each of the intervals [c,_, , c,] leads to 

ab n ack n 
| p(t) dt = >> | p(t) dt > >> 17/d, , 
Fa k=l Jeg-; k=1 
where d, = ¢, — C:-; > O. Since } a we d, = b — a, this implies (2) for an arbitrary n 


if it is ascertained that 


IV 


n'/ > d, >. the, . (4) 
k=1 k=1 
But (4) is true, since it merely expresses the fact that the harmonic mean of n positive 
numbers d, cannot exceed their arithmetical mean.’ 
Accordingly, it is sufficient to prove (2) under the following two (simultaneous) 
assumptions: n = 1, and [a, b] is a complete primitive interval of some solution x(t) of 
(1). Then, if a = 0 without loss of generality, the assertion (2) reduces to 


3 


b | p(t) dt > 1, (5) 


and the assumption of (5) is the existence of a solution z(t) for which both x(t) # 0 
and 2’(t) + 0 hold on the open interval 0 < ¢ < b and etther x(0) = O and z’(b) = 0 
or x'(0) = 0 and x(b) = O. But it will be clear (by interchanging “past” and “‘future’’) 
that it will be sufficient to prove (5) for the first of these two alternative cases. 

In addition, since x(t) ¥ 0 for 0 < t < b, and since the solution x(t) can be replaced 
by — z(t), it can be assumed that z(t) > 0 for 0 < ¢ < b, while z(0) = 0. Since, as 
pointed out above, the graph of x = z(t) always turns its concavity toward the f-axis, 
and since z’(t) 0 for 0 S t < b [while 2’(b) = OJ, it follows that, if 0 < t < b, the 


See, for instance, G. Pélya and G. Szegi, Aufgaben und Lehrsdtze aus der Analysis, vol. 1, 1924, p. 50. 
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ordinate of the graph of x = x(t) is positive and increasing, while its slope is positive 
and non-increasing. Since the slope at ¢ = 0 is 2’(0) > 0, and since z’(b) = 0 prevents 
that x(t) be linear on the whole of [0, 6], it now follows from x(0) = 0 that x(b) < bx’(0). 

On the other hand, since x«'(b) = 0, integration of (1) between ¢ = 0 and ¢ = b shows 
that 


ab 
x0) = | pal at. (6) 
But since x(t) increases from x(0) = 0 to x(b) > 0 on the f-range of (6), and since p(t) = 0, 
the representation (6) of 2’(0) implies that, unless p(t) vanishes identically, 


ab 
0 < 2x0) < 2(b) | p(t) dt. 
0 


Hence (5) follows from the inequality x(b) < bx’(0), found at the end of the preceding 
paragraph 
4. This completes the proof of the italicized assertion. The following fact is a corollary: 
If (1), where p(t) = O anda St S b, possesses some solution for which x(t) and x(t) 
together do not have more than one zero on the open interval a < t < b, then 


ab 


| p(t) dt > 4/(b — a). (7) 

In fact, since (7) is the case n = 2 of (2), it is sufficient to ascertain that the assump- 
tions of (2) are satisfied by n = 1 or n = 2 according as x(t) and x’(t) together have 
no or a single zero on the interval a < t < b. But this is clear for reasons of concavity 
(see above), since the whole of [a, }] is a primitive interval of «(¢) in the first case, while 
[a, b] consists of two primitive intervals of x(¢) in the second case. 

5. It is well-known that a result of Liapounoff’ is substantially equivalent to the 
following assertion: If p(t) = 0, then (7) must hold whenever (1) has on [a, b] a solution 


x(t) satisfying 
ria) = 0. x(b) = 0, x(t) > 0 for a<t< b. ©) 


[t is also known’ that, under the assumption (8), the 4 is the best absolute constant in 
7). A simple proof of the fact that (7) is necessary for (8) was given by Shukovski (a 
contemporary of Liapounoff)* and recently by Borg.* 

Clearly, this result is a particular case of the italicized corollary (n = 2) of the 
general inequality (2). In fact, the assumption (8) means that the graph of x = 2x(f) 
is a convex arch over [a, b], and so z’(¢) has just one zero on [a, b]. But (7) turns out to 
hold also for boundary conditions more general than those of (8). For instance, (7) 
must hold also if (1) possesses a solution x(t) which satisfies the boundary condition 
x'(a) = 0, x’(b) = O but has only one zero on [a, }]. 

The general result, expressed by (2) for an arbitrary n, suggested itself by Shukovski’s 


Actually, Liapounoff considered only the case of an equation (1) having a periodic coefficient 
p of period b — a; see A. Liapounoff, Comptes Rendus 123, 1248-1252 (1896). 

se. R. van Kampen and A. Wintner, Amer. J. Math. 59, 270-274 (1937). 

‘See J. L. Geronimus, Alerander Michailowitsch Ljapunow, Berlin, 1954, pp. 52-53. 

5G. Borg, Amer. J. Math. 71, 67-70 (1949). 
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and Borg’s proof of (7) for the case (8). In fact, even the use of (4) is” between the lines 
of that step in Borg’s proof of (7) in Liapounoff’s case which leads from [a, c] and [c, }] 
to [a, b], where c is defined by x(c) = max x(t). 


MAXIMUM SPEED IN STEADY SUBSONIC FLOWS* 
By CHIA-SHUN YIH (University of Michigan) 


The purpose of this note is to show that the maximum speed in a singularity-free 
region of any steady subsonic flow must occur on its boundary, provided the flow is 
isentropic and irrotational. The corresponding result for irrotational flows of an in- 
compressible fluid is well known. 

Since the flow under consideration is irrotational, the square of the speed g is given by 


q = ¥.%,; (1) 


in which ¢ is the velocity potential, the summation convention is used, and, with Cartesian 


coordinates, 


is the jth component of the velocity vector—with the usual minus sign suppressed for 
convenience. The equation of continuity is, for steady flow, 


in which p is the density. The equations of motion are 


‘ 


eee 
(In p).; = a2 (G + 2923); (3) 
if x3 is measured in a direction opposite to that of the gravitational acceleration g, and 
c is the local speed of sound. Finally, the Bernoulli equation can be written in the form 


(with p indicating pressure) 


j=, - a3 (q° + 2923) (4) 

Po 2 
in which k is the ratio of the specific heat at constant pressure to that at constant volume, 
and the subscripts zero indicate that the quantities involved are taken at some point 
of reference, where the gas is at rest and from which x, is measured. 

For a proof, it is sufficient to show that the Laplacian of g° is non-negative. Whereas 
this is easily done for irrotational flows of an incompressible fluid, it cannot be readily 
established for the flows under discussion. Professors D. Gilbarg (through his publi- 
cations) and C. Truesdell (by oral communication) have indicated to the writer that 





6G. Borg, Amer. J. Math. 71, 68 (1949); see also P. Hartman and A. Wintner, ibid., p. 209. 
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the non-negativeness of Vg" can be replaced by that of 
(bi i(Q) i: 
if the matrix 6,; is positive definite. It will be seen that their suggestion provides one 


of the key steps to the following proof. 
Taking the Laplacian of Eq. (1), one has 


(T).is = @.i.i) 8s = Weiss H Wiwiss - (5) 
But Eq. (2) can be written as 
¢.i:i + (In p).e,; = 90 
so that 
Qi Qi —[(In p) ee. = —(n p).¢.i¢.; — (In p)i¢.¢.; 

and Eq. (5) becomes 

(Qi = 2.49.4; — (In p).¢.:¢,.; — 2(In p)i:9..9.; - (6) 
3y virtue of Eqs. (3) and (4), one has 


l > k — ] 2 > Cod 
-2(In p) .i¢.¢ 2D)... F Ga GH ges) (G+ 2945) .10.08.; - (7) 


2¢.¢.5: = (7). #90 
then g cannot be a maximum. Therefore, one can concentrate on the case 
29.0.4: = (¢)., = 0. (8) 


In this case Iq. (6) becomes, after substitution from Eqs. (7) and (8), 


I , 29°(k — 1) 
(VQ). — 3.0440) 65 = 20,110.55 + — 53559, ; 
ol 
2g°(k — 1) ses 
b.,(4q 2¢ 4¢.5; += a - 6,36;.9..¢.; = positive (9) 
in which 
. l - i 
b, 6; —-3¢.¢.;, and 6,;,; = the Kronecker delta. (10) 
2 
Now since the flow is subsonic, 
CoPia o | 
c iad 
Then obviously if 
bi, = ="* 3,, -sa¢e.¢ (11) 
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is positive definite, b;; must also be. The proof will therefore be complete if one can 
demonstrate that b/; is positive definite, i.e., the bilinear form b/; u;u; is positive definite— 
u, , and wu, being any three real numbers. But 


Uy 


b! uu; = [(wie.2 — Ue.) + (ee.3 — Use.2)) + (Us¢.1 — WMy.3)° (12) 
Thus b!; and therefore b;; are positive definite, and Eq. (9) shows that the maximum 
value of g must occur on the boundary of any singularity-free region. 

Since the Bernoulli equation is 
I; ) ( 
I + A 4 gx; = constant 


y-15* ¥" 


and the equation for isentropic change is 
-k ~k 
Pp = PoPo = constant, 


minimum pressure, minimum density, and (from the equation of state for perfect gases) 
minimum temperature are always associated with the maximum value of 
q + 2923. 

If now one takes the Laplacian of g° + 2g2; , one can show by a procedure strictly similar 
to the foregoing one that the maximum value of q’ + 2g7; must occur on the boundary 
of any singularity-free region. Thus for any such region of a steady subsonic flow, the 
minimum density, pressure, or temperature must occur on the boundary, provided the 
flow is irrotational and isentropic. 

The boundary of a singularity-free region is, of course, not necessarily a_ solid 
boundary. However, in the case of an ambiently uniform flow past closed solid bodies 
placed in an infinite fluid, one may compare a streamline at infinity with another not 
entirely at infinity. Both streamlines can be considered to originate and terminate at 
the same potential planes (on each of which ¢ is constant). Since there is variation of 
speed on the latter streamline, the maximum speed thereon must exceed the uniform 
fluid speed at infinity. This fact, plus the theorem for the speed just proved, shows 
that, for any ambiently uniform flow of the kind under discussion, the maximum speed 
must occur on one of the solid boundaries present in the flow if it is free from singular- 
ities—since the entire fluid space bounded by infinity and the solid boundaries can be 
taken to be the region under consideration. 

A flow with uniform ambient velocity U past a solid body is physically the same 
as one caused by that body moving with velocity — U in an otherwise quiescent infinite 
fluid, i.e., the distributions of pressure and density relative to the body must be the 
same for both cases. Thus, if gravity is ignored, the minimum density, pressure or 
temperature in a subsonic flow caused by a body moving steadily in an otherwise quiescent 
fluid must occur on the boundary of that body if the flow is isentropic, irrotational, and 


free from singularities 
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A NOTE ON THE NUMERICAL INTEGRATION OF A CLASS 
OF NON-LINEAR HYPERBOLIC EQUATIONS* 


By R. BELLMAN (The RAND Corp., Santa Monica, Calif.) 
anp I. CHERRY anv G, M. WING (Los Alamos Scientific Laboratory, University of 
California, Los Alamos, N. M.) 


1. Introduction. Since the equations of hydrodynamics are non-linear, a computa- 
tional solution of these equations usually involves a numerical integration. This numerical 
integration is particularly difficult to perform in the presence of discontinuities of the 
solution called ‘‘shocks.’’ Although there are a number of devices available for obtaining 
the solution in the neighborhood of the shock, notably by von Neumann, Richtmeyer 
and Lax, these have the unfortunate property of altering the behavior of the solution 
in other regions. If conventional finite difference schemes are employed, instabilities 
arise in a number of cases. 

In this paper we propose a new method of numerical integration which seems well 
suited to non-linear hyperbolic equations of the form 


u, = gluju, + h(u, x, 0, u(x, 0) = f(z). (1) 
To illustrate the method, we consider the pseudo-hydrodynamic equation 
u, = —Uu, , u(x, 0) = f(x), (2) 


which exhibits the shock phenomenon and which possesses the great merit of an analytic 
solution. We can thus compare the analytic solution with that obtained from the appli- 


cation of our numerical technique. 


2. The equation u, = —uu,. The equation appearing in (1.2) has the solution 
u(x, t) = fix — u(x, Ot). (1) 
Since 
'(E) . ; 
= soon hats E=2-— u(z, bt, (2) 


1+ ¢f(E£)’ 
it is clear that u, becomes infinite whenever 1 + ¢f’(£) = 0. This is the shock phe- 
nomenon. Numerical integration is difficult in the neighborhood of this shock. 
In place of the usual finite difference method based upon a fixed (x, t)—grid, consider 


the following scheme: 
u(x, t+ A) = ulx — u(z, OA, ft], u(x, 0) = f(x), (3) 


where ¢ assumes the values 0, A, 2A, --- . Expansion of the two sides in powers of A 
reveals that (3) is equivalent to (1.2) to first order terms in A. 

The method does not involve the strict notion of an (2, ¢) grid, since the x — u(r, A 
values will not be distributed in any regular fashion. In integrating the equation, a fixed 
z-grid is used, x = --- , —26, —6, 0, 6, 26, --- . When the values u(né, ¢t) are known 








*Received April 23, 1957; revised manuscript received June 18, 1957. Work done under the auspices 


of the U. S. Atomic Energy Commission. 








182 NOTES [Vol. XVI, No. 2 


for fixed ¢ the quantities 
u[néd — u(nd, OA, t] = u(nd, t + A) 


are found by interpolation in the x direction. There is no constraint on the relation 
between 6 and A, and 6 may be chosen as small as desired, subject to space limitations 
on the digital computer and the cost of computing time. 

One important feature of the recurrence scheme in (3) above is that uniform bounded- 
ness of u(x, t) over x fort = kA is automatically passed on to u(x, t) att = (k + 1)A. 
3. An example. We have integrated (2.3) using for the starting condition 


(0 , -10<2 <5 -05; 
9/12(2 + 1/2), —05 Sz Ss —0.26; 
f(a 0.1 , —-0.26 sx s 0.26; (1) 
5/12(1/2 — zx), 0:25 s 2s 035: 
0 : 07 S24 Ss 10: 
f(x) f(x + 2), —o <7 < o@ 
We have used six-point Lagrangian interpolation, x points at intervals of 6 = 0.02, 


and. various values of A. Some results are shown in the table. It should be noted that 


even with relatively large A the discontinuity in u that takes place at t = 2.4 is handled 
quite well. Other starting conditions have also proved successful. 

4. Discussion. Since the method consists of integrating along the characteristic, 
good results are to be expected. However, we have also integrated equations of the form 


u, = gu, + hu, 2, (1) 
u(x, O) = f(r), —o2o < ££ < @, 
using the difference scheme 
u(x,t + A) = u(x — glu(x, t)]A, 0 
+ hlu(x, t), x, tA, (2) 
uz.) = f(x), =—oO< Fs <x wo, 
where the characteristic curves are not simply x — g(u)t = k, and have obtained ex- 
cellent results. 
5. Extension. The method is applicable to the integration of the system 
u, G(u,v, x, yu, + Hu, v, x, yu, + glu,v, x, y), (1) 
v, M(u,v, x, yu, + N(u,v, x, yu, + h(u,v, 2, y). 
The corresponding difference scheme is 
u(x, y, t+ A) = usa + GA, y + HA, t} + gA (2) 
v(x, y, t+ A) = via + MA,y+ NA, t} + ha. 


It is readily verified that (5.2) is equivalent to (5.1) to terms of order A. The inte- 
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gration proceeds as in the previous case, a grid in the (2, y) plane being chosen and 
interpolation among these grid points being used to advance the integration in time. 
However, attempts at extending the integration method in an obvious way to general 
systems of hyperbolic equations in one space variable have not been successful. While 
the method presented here is marginally stable, the attempted generalizations have 
suffered from instabilities. 
TABLE 1 


Values of u(z,t) att = 2.4 


x Exact u Computed u 
A = 0.02 A = 0.10 A = 0.20 A = 0.30 
— 48 | 0.004167 0.003713 0.003657 0.003584 | 0.003508 
— 0.0208 0.020805 0.02056 0.02024 0.01991 
— 32 0.0375 0.03738 0.03694 0.03637 0.03576 
— oo 0.0542 0.05399 0.05336 0.05253 0.05165 
<2 0.0708 0.07063 0.06981 0.06869 0.06745 
— .08 0.0875 0.08748 0.08621 0.08484 0.08255 
00 0.1000 0. 10022 0. 10034 0.10000 0.10072 
+.08 0.1000 0.09994 0. 10000 0. 10000 0. 10002 
+.16 0.1000 0. 10000 0. 10000 0. 10000 0. 10000 
24 0.1000 0. 10000 0. 10000 0. 10000 0. 10000 
+ .32 0.1000 0. 10002 0. 10000 0. 10000 0. 10000 
+ 40 0. 1000 0. 10028 0. 10008 0. 10000 0. 10009 
+ 44 0.1000 0.09924 0.09952 0. 10000 0.10111 
+ .46 0.1000 0.09785 0.10351 0.09884 0.08966 
48 0.1000 0. 12222 0.08834 0.06096 0.04563 
50 0.00000 0.00000 0.00000 0.00000 
+ .52 0.0000 0.00000 0.00000 | 0.00000 / 0.00000 
+ 54 0.0000 0.00000 6.00000 0.00000 0.00000 
56 0.0000 0.00000 0.00000 0.00000 0.00000 


ON THE EQUIDISTRIBUTION OF PSEUDO-RANDOM NUMBERS* 
By J. N. FRANKLIN (California Institute of Technology) 


1. Introduction. In the paper [1] J. Moshman discusses a method for generating 
pseudo-random numbers. As Moshman points out, his method is an adaptation to decimal 
computers of a process used on a binary computer by J. Todd and O. Taussky Todd [2]. 
Another adaptation to decimal computers was made by M. Juncosa [3]. The method is 
used widely for computations in nuclear physics. For example, the method is used by 
E. Leshan [4] to study the transport of neutrons. 

In this method one chooses an initial number x) and a multiplier N. Each of these 
numbers is represented in a digital computer by a word with a fixed number of digits, 
say s digits. The product Nz, contains 2s digits, of which the s least significant digits 
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are used to represent x, . Now 2, is represented by the least significant s digits of the 
2s-digit product, Nz, , and so on. 

In a digital computer the number 2, is rational, and the sequence 2 , 2 , % , --~ is 
periodic. In the present paper the sequences 2 , 2; , %2 , «++ are considered in which 2 
is zrrational. By means of the ergodic theorem it is shown that almost all of these sequences 
are equidistributed. 

2. Application of the ergodic theorem. Let us normalize the number x, by requiring 


0. ey < f, (1) 

We set N equal to an integer > 1. Then 
z, = Nrz,-, — [Nz,-:] va = 1,2, --:-), (2) 
where [Nx] stands for the greatest integer < Nz. In this notation the Todds, Juncosa, 


and Moshman choose, respectively, 


vy = 27° N=5 
t%, = 107’, N = 3°"! 
m4k+l 


Lo = | eee N=7% ; 


where s is the number of digits in a word of the computer, and k is some positive integer. 
One may consider more general sequences by fixing a number @ and setting 


x, = Na,-, + 6—-— [Na,-1 + 9] (mn = 1,2, ---). (3) 

Without loss of generality, we suppose 0 < 6 < 1. When @ = 0, we have the definition (2). 

The sequence 2) , x; , «+ is said to be equidistributed if, for every fixed a, b in the 

rangeQO<a<6b< 1] 

l 

- a l~ b—a as n—-o. (4) 

When .V = 1 in (3), we have the case considered in the classical paper [5] by H. Weyl, 

who showed that for every value of x, the sequence (3) where N = 1, is equidistributed 


if and only if @ is irrational. 

By Weyl’s result the sequences (3) where N = 1, are equidistributed for all x) and 
almost all 6. If the integer V is greater than 1, we shall show that the sequences (3) are 
equidistributed for all @ and almost all x 

Theorem. Let f(x) be an arbitrary integrable function in the sense of Lebesgue. Let 


0 < x, < 1, and let a sequence 2 ,z, , 2 , --- be formed according to (3), where N is a 
fixed integer > 1 and @ is a fixed number with 0 < 6 < 1. Then for almost all x 
~~ . F ‘ 
pm f(a) | f(x)dr as no, (5) 
N k=0 Ji Z 


The equidistribution (4) for almost all xz) follows at once from the theorem if f(x) 
is defined to be the characteristic function of the interval a < x < b, which takes the 
value 1 or O according as = lies within or outside of the interval. 

In order to prove the theorem we make use of a statement of the ergodic theorem 


due to F. Riesz [6]: 
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“Tet a measurable set 2 be given, of finite or infinite measure, the corresponding 
measure and integral being defined according to Lebesgue, or more generally, by means 
of a distribution of positive masses. That being the case, let us designate by 7 a point- 
transformation which is single-valued (but not necessarily one-to-one) from 2 onto 
itself; and let us suppose that 7’ conserves measure in the sense that, F being a measur- 
able set, 7'Z its transform, and E’ the set of points P whose images appear in TE, the 
sets E’ and TE have the same measure. Then, if f,(P) is an integrable function and 
f.(P) = f,(T*""P), the arithmetic mean of the functions f, , f. , --- , f, converges almost 
everywhere, as nm — ©, to an integrable function g(P) which is invariant (almost every- 
where) under 7’.”’ 

Riesz gives an elegant proof of this form of the ergodic theorem, which ‘s more general 
than G. Birkhoff’s original form [7] in that the transformation 7’ is not required to have 
a unique inverse. After completing the proof, Riesz writes without elaboration: 

“Let us add finally that in the case where Q is of finite measure, it follows by inte- 
gration term by term (which is permitted in this case because of the uniform integrability 


of the terms) 


/ oP) = [ £P).” 


A proof of this statement is written out at the end of this paper. 
In our case the set 2 consists of the interval 0 < 2 < 1. The transformation T is 


Tx = Nx + 6 — [Nx + 8]. (6) 


Each point x in Q is the image under 7’ of exactly N points: 


pkta-60 (k=0,---,N—-1 if «>@, ‘ 
1 a ae (7) 
’ lk =1,---,N Zz 2< @ 
To show that 7 preserves measure, it is sufficient to show that 
1 1 
[ f(Tx) dx = [ f(x) dx (8) 
“0 “ O 


for all integrable functions f(x). In fact, if f(a) is the characteristic function of a measur- 
able set TE, then f(7'x) is the characteristic function of the set E’ of all points whose 
images lie in 7'E, and Eq. (8) states that the sets 7’F and E’ have the same measure. 
By the periodicity of the function x — [x], we have 


al 


[ f(T) ax = [ "(Na + 6 — [Nx + 6]) de 


Je 
N 


1 a. 

V | fly + 6 — [y + 6]) dy 
= [ fy + 6 — [y + 4]) dy 
al+6 


= l fly — ly}) dy 


1 
= [ fly) dy, 


which establishes the required formula (8). Now we can prove our theorem. 
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Proof of the theorem. By the definitions (3) and (6), 
= Tx, (k =0,1,2, ---). (9) 


By the ergodic theorem, there exists an integrable function g(x.) defined almost every- 
where for 0 < 2) < 1 such that 


1 n-1 
- 7 f(a.) ~ (a) ae.as n> @, (10) 
¢(Txo) = g(t) 2a.€., (11) 


where “‘a.e.”’ means ‘“‘almost everywhere.’”’ To prove the theorem we must show that 
the function g(x) is almost everywhere equal to a constant, namely the average value 
of f: 

al 


g(t) = f(x) dx a.e. (12) 
Let the integrable function ¢(t) have the Fourier series 
>> c, exp (2kmit). (13) 
k=-@ 


According to the theorem of Fejér and Lebesgue [8], p. 415, this series is summable by 
arithmetic means to the sum ¢(t) for almost all ¢. 
We shall show that 


c, = exp (—2kri6)ey, (k=0,+1,+2,--:). (14) 
By the invariance (11) of g, we may write 
On = | g(Nt + 6 — [Nt + 6]) exp (—2Nkrit) dt 
=x | o(t + 6 — [¢ + 6) exp (—2kwit) dt. (15) 
By the periodicity of the exponential function, 
= | g(t + 6 — [t + 86]) exp (—2kmit) dt 
al 
= | g(t + 6 — [t + 6]) exp (—2kwit) dt 
al 
= | g(t) exp [—2kri(t — 6)] dt = exp (2kwi8)c, . (16) 
This establishes (14). 
By iterating (14) r times and taking absolute values, we find 
ce, | = |c,. |, with m= N’k (ry = 1, 2, ---). (17) 


But the Fourier coefficient c,, ~ 0 as m — + o. Keeping k fixed in (17) and letting 
r— o, we find 
c¢=0 if k x0. (18) 
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Therefore, the Fourier series (13) contains only the constant term c) . Then the sum by 
arithmetic means equals ¢c, , and by the Fejér-Lebesgue theorem 


Af) =c a.e. (19) 


It remains only to show that the constant c, has the value 


1 
c= | Hla) ae. (20) 
This follows immediately from the statement fg(P) = Jf,(P), which is appended to the 
ergodic theorem. In our case, by (19), 


1 


1 

Co = / g(t) dt = [ f(x) dx. (21) 
0 “0 

This completes the proof of the theorem. 

3. Proof that fy(P) = Jf,(P). As Riesz suggests, this statement, which is appended 
to his generalized ergodic theorem, may be proved by Lebesgue’s convergence-theorem 
if the underlying space Q has finite measure. 

For any number \ > 0 we may define a bounded function g (P) equal to f,(P) or 
to 0 according as | f,(P) | < Aor) f,(P) | > ’. We may then define the remainder h,(P) = 


f,(P) — g,(P). Let an arbitrary number ¢ > 0 be given. Since f,(P) is integrable, we may 
choose a number A = XA(e) so large that the corresponding function h,(P) satisfies the 


inequality 


a 


| | h(P) | <e. (22) 


72 


We now apply the ergodic theorem separately to the functions g,(P) and h,(P), 
obtaining limits 


lim ; > g(P) =x(P) ae, (23) 
ner kel 
lim * > h(P) = WP) a.e. (24) 
sill k=1 


Addition of these equations gives 
g(P) = x(P) + WP) ae. (25) 


Since every mean value appearing on the left-hand side of (23) has absolute value 
< ) for all P, and since m(Q) < ©, we may apply Lebesgue’s theorem of bounded 


convergence to obtain 


- 1 : io >) Oye 
lim Df olP) = [ x). (26) 

Since fg,(P) = Jg,(P) for all k, (26) gives 
[ oe) = [ x. (27) 


Even though the function h,(P) is integrable, it is not obvious that all of the mean 
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values 


are majorized by a single integrable function H(P). Therefore, to formula (24) we apply, 
not Lebesgue’s theorem, but Fatou’s theorem. Using the inequality (22), we find 


. ah n 
| ¥(P) | < lim | F >> h.(P) 
a t l 


“ 
nc 


lA 


: F _. my > 
- . ) 2 | h(P) | = | h(P)| <e. (98) 


Combining our results, we find 


“ 


[ {(P) - [ ¢(P) | | gi(P) — [ x(P) + [ h(P) — [ ¥(P) 


« 


2 


</ | a) - | x(P) | + | nae? te [ | WP) 


< 


(29) 


<O+e+e. 


Letting « — 0, we obtain the required result {f,(P) = fe(P). 
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SIMULTANEOUS INVARIANCE OF GENERALIZED SPHERICAL 
HARMONICS UNDER THE OPERATIONS OF TWO ROTATION 
GROUPS* 


By R. N. DDHEEDENE (Cornell University) 

Abstract. A method is found for evaluating the coefficients of a sum of generalized 
spherical harmonics so that the sum will be simultaneously invariant under two rotation 
groups. The coefficients for a sum of ordinary spherical harmonics invariant under 
each individual group must first be known. 


*Received May 15, 1957. 
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If a rigid molecule is situated in a lattice, free to rotate but not displace, the mole- 
cule’s potential energy V, a function of its rotational position, must have certain sym- 
metries, dictated by the symmetry of the molecule and the symmetry of the lattice. In 
this paper a consequence of these symmetry properties of V is discussed. 

Let the position of the molecule be described by the three Eulerian angles (6, ¢, y), 
of which @ and ¢ are the polar coordinates of a z-axis fixed in the molecule with respect 
to an (X, Y, Z) system of coordinates fixed in the lattice, and y describes the rotation 
of the molecule about z'. Let there be a group G of rotations of (X, Y, Z) describing the 
symmetries of the lattice, and a group H of rotations of (x, y, z) (the molecule-fixed 
coordinates) describing the symmetry of the molecule. Suppose V (6, ¢, ¥), the potential 
energy, is expressed as a series of generalized spherical harmonics (often called the 
symmetric top eigenfunctions), 

» J J 
V0eW=> dX} YD Ch”yT (6, ¢, v). (1) 
J=0 K=-J M=-J 
Then the requirement that V be invariant under each operation in G or H will place 
certain restrictions on the constants C/:” (the higher the symmetry expressed by G 
and H, the more of them will be zero). We wish to determine these restrictions. 
Instead of using the real triplet (@, ¢, y) to represent Eulerian angles, we will use 


the complex doublet (or spinor) (£, 7) where 


6 
~ = exp [—i(¢@ + ¥)/2] cos5 , 


7 (2) 
oie, & 
n = exp [7(@ — y) 2]sins , 
or even more simply, the matrix, r.. where 
t* - 
r=({°* 7”). (3) 
-—9 § 

In this notation, if (x, y, z) is given by T, and (2’, y’, 2’) by T. , each with respect to 
(X, Y, Z), and (x’, y’, 2’) by T with respect to (x, y, z), then the relation between the 


angles is 


r,=TP,F, (4) 


an easy result, but one involving some arithmetic’. 

We will write the generalized spherical harmonics as Yf'” (T°) where T indicates 
the argument (6, ¢, ¥) of Y. 

Note that we may consider (4) as giving the transformation under I of the Eulerian 
angles T, into f. . That is, I represents a rotation R (1) in H, under which Tr, — I, 
gives the transformation of the Eulerian angles describing (x, y, z) with respect to 
(X, Y, Z). An important result is that under such a transformation, 


J 
RIV)YsS"(r,) = VF") = Do Yer) ¥7'"(r) (5) 
S=-J 


ISee Goldstein, H., Classical Mechanics, Addison-Wesley, 1953, p. 107 for an excellent figure show- 


ing the Eulerian angles. 
2Takahashi, T., J. Phys. Soc. Japan, 7, 307 (1952). 
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gives the transformation of the generalized spherical harmonics’. Letting Y,(T) be a 
° ° 7K. Miv an ° 
square matrix with elements Y#'“(T), (5) may be written as 


R(T)YAT,) = Y,(T.) = Y,(T,) Y,(P). (6) 


The transformation matrix elements for ordinary spherical harmonics are the same as 
those for generalized spherical harmonics; using essentially the same notation, 


R(T) Y5(8; 561) = Y5(O. .¢2) = LL Y5(@ , ob) ¥F"*(P) (7) 


is a well known result’. 

The above equations express transformations under a rotation R(T) in H. For a 
transformation R(T’) in G, the situation is slightly different. If we let R(T’) (X, Y, Z) = 
(X’, Y’, Z’), and T, , T, locate (2, y, z) with respect to (YX, Y, Z), (X’, Y’, Z’), respect- 
ively, then as we see by applying (4) and considering (X’, Y’, Z’) as the (X, Y, Z) of 
equation (4), we have 


r, = (r’)'T, (8) 
(exponential-1 indicates inversion) replacing (4), so that 
R(T’) Yi(T:) = YT.) = Yi((1")’) YP) (9) 
gives the transformation of the generalized spherical harmonics for elements R(T’) of G. 
Now let the sum in (1) be invariant under each R(T;) in H and each R(T’) in G. 
Then, by (6) and (9), we require for each R(T;), 


Zz. CF gt fame | it Ri r,) >. c* ef af r) 


oar J.K.M (10) 


J,K,M,8 
and, using the fact that if R[(I’)~'] is in G, then, since G is a group, R(T’) is in G, for 
each R(T‘), 
es ee RD > SI 
J.K,M JK,M (11) 
= :% cr uf y; M r) y* S(T’), 
J.K.M.S8S 


Equating the coefficients of Yf'“(T) in (10) and (11) gives 


Cr = VeCre ye (ry) = Verh YF (Pr) (12) 


for each I, in H, T inG. 
Let constant vectors h; , g, be known, such that 


YT hy = h, ’ VM a = (13) 


for each T,; , Tj . It is reasonable to assume such vectors known, since, using group 
theory to determine the number of identical (A,) irreducible representation occurring 


‘Wigner, E. P., Gruppentheorie und thre Anwendung auf die Quantenmechanik der Atomspektren, 


Edwards Brothers, 1944. 
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in the 2J + 1 dimensional representation of the group, one can, at least in principle, 
find coefficients h* so that the sum 
J 
DY AF Y5(6, 4) (14) 
Ke-J 


is invariant under a group H (the same comments obviously hold for g, G)*. And, by 
(7), this gives, for each T; in H, 


Pi 
> AT Y5(0,¢) = R(T.) > kT Y5(8, W) = 
Ku-J K (15) 
= > AsY5(0,¢)¥;'*(r,). 
K.S 


Equating the coefficients of Y7 (6, @) in (15) gives (13). 
We now have the solution to our problem; setting 


CP = hy'gy* (16) 
makes (1) invariant under both G and //, as we see by showing directly that (16) satisfies 
(12), using (13) and the fact that Y, is Hermitean: 


DEP EVE) = Lehrer VP) 


— Avgk« vK,M 


Q = ° —_ 
x Qs J (1 7) 
> hs or* ¥3'*(r,) 


WS rS.K Ta 
7c "Y; (T’;) 


S 


KypM _ ~vK.M 
gy hy = C; ° 


Generalization of this theorem to finding appropriate C*'” so that the sum (1), 
rather than being invariant, will have specified transformation properties, i.e. belong 
to a representation other than A, , is obvious. 

As an example the first non-zero Cf'™ will be found for the case of a tetrahedron in 
a cubic lattice (e.g. solid methane). For the cubic group, A, occurs once in the one 
(J = 0) dimensional representation, not at all in the representations corresponding to 
J = 1, 2, 3, and once in the nine (J. = 4) dimensional representation. Similarly, for the 


tetrahedral group, A, occurs once in the representations corresponding to each of J = 0, 


3, 4. It is trivial that Cy # 0. Since no A, representation occurs for J = 1, 2 in either 
group, C,; = C, = 0. Since g, = 0 (corresponding to the cubic group), C; = 0 although 
h, ~ 0. Now the sum 

fm \1/2 ~ \1/2 

(=) P? (cos 6) + (3) (2)'”’P{ (cos 0) cos 4¢ (18) 


of fourth order ordinary spherical harmonics is invariant under the cubic group‘ and 
therefore under the lower-svmmetry tetrahedral group. Thus 


5 1/2 7 1/2 5 1/2 |7 
a = h, = (3) ’ 0, 0, 0, (5) ’ 0, 0, 0, (3) (19) 


‘Bethe, H. A., Ann. der Physik, ser. 5, 3, 133 (1929). 
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and (16) gives 


5 5 os ; — so 
py (Ya + er +4 . + Fe *) 
2 (20) 
ive ' } 0 } ‘ ' rors ane re 
— (35 4 T j + 1 be 4 : + 12 4 
as a properly invariant function. 
The author wishes to extend his appreciation to Professor R. Bersohn, without 
whose assistance this work would have been impossible. 


ON THE MOTION OF A SIMPLE PENDULUM* 


By BORIS GARFINKEL (Ballistic Research Laboratories, Aberdeen Proving Ground, Md. 


Abstract. The vanishing of the tension in a simple pendulum supported by a flexi- 
ble cord causes the particle to pass from the circular to a parabolic trajectory. The 
number and the nature of such transitions are related here to the value of the initial 
energy. 

1. When the initial energy of a simple pendulum lies in a certain interval, the tension 
vanishes at some instant of the motion. Then, if the support is provided by a flexible 
cord, the particle passes from the circular to a parabolic trajectory. The number and 
the nature of such transitions are shown here to be precisely related to a dimensionless 
energy parameter, ~. Despite the intrinsic interest and the relative simplicity of this 
motion, it does not appear to have been treated in the literature. 

Let | be the length, m the mass of the pendulum, 7 its distance from the point of 
support, @ the angular coordinate measured from the downward-drawn vertical line, 
and g the acceleration of gravity. The constraint 


can be replaced by the condition 
Al — r) = O,7 (2) 
where ) is a multiplier vanishing if 1 > r and admitting a non-zero value if 1] = r. The 
Lagrangian of the system, 
L = 3m(r? + 7°60") + mgr cos 6+ Xl — 7), (3) 
leads to the differential equations of motion, 
mr’ — r0? — g cos —) +X = 0, r0° + 2r6 + gsin 6 = 0, (4) 


which together with (2) and (1) determine the functions r(¢), 6(¢), \(t) when the initial 
conditions are prescribed. From (4.1) the multiplier \(¢) can be identified with the 


*Received July 22, 1957. 
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tension of the cord. Let the initial conditions now be: 
(0) = 1, 6(0) = 0 
} ’ ( J ’ (5) 
r'(0) = 0, 0°(0) = (2E/ml?’)'”, 


E being the initial energy measured from the lowest point. Then the particle moves 


in a circular are r(t) = 1 until \ vanishes. The integrals of the motion, deduced from 
(4), are 
(0°/w)* — 2 cos 6 = 2(2c° — 1) = 3é, sin (6/2) = min (1, ¢) sn [wt/max (1,c)], (6) 


where the constants w, c, are defined by 
w = g/l, c = E/2 mg l, & = 2(2¢ — 1), (7) 
and the modulus of the elliptic function sn is 
k = min (c, 1/c). (8) 


The dimensionless parameter £ is seen from (7) to be proportional to the energy measured 
from the horizontal configuration r = 1, 6 = w/2; its range is — 2/3 < — < @. 

2. At the instant ¢ = ¢, of the vanishing of the tension the coordinates and their 
derivatives, obtained from (4.1) and (6.1) with the substitutions \ = 0, r = 1, rr = 0, 
are given by 

cos 6, = —€& n=l 
i $) 1 ’ (9) 

6;/o = ¢'”, r, = 0. 

Hence @, and 6; are real if and only if & lies in the range 
O0<i<l, (10) 
and 7/2 < 6, < r,0 < 6; < w. Att = ¢, there occurs a transition from the circular to 
a parabolic trajectory. The latter corresponds to the solution of (4) with A = 0, r < J, 

and the initial conditions (9), and is represented by the equations 


(r/l)? = 1 — (wr)* (EC — &)]'? — (r)*/4, (11) 
rsin 0/1 = (1 — #)'? — #°@7), 
where 
t= l = t; . (12) 


When the particle re-enters the circle, r again assumes the value r = 1, and (11) yields 


2 


wr = 4[e(1 — &)]'”, 
9, = Ir — 36, , = 
i] vs 30 lr. l (13) 
6; = 6)(—3 + 122° — 8#'), 

r; = Slw[é(1 — &)]”. 


Such an alternation of the trajectory between the circle and a parabola we shall, for the 
sake of brevity, designate by the term “flip.” At the instant t = ¢, of re-entry of the 
circle the ideal cord, assumed to be inextensible and infinitely strong, enforces the 
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constraint (1) by generating an impulse equal and opposite to the radial momentum 

m’. While the latter is being annihilated, the angular momentum mr°@ is conserved, 

and the energy is thus diminished by the quantity mr™*/2. The corresponding diminution 
-/ 


of — can be calculated from (13.4) and the definitions (7); the new value ¢’ can then be 
written as a function, f(£), in the form 


él —£)])*, for O<é <1. (14) 


Of course, if — lies outside the range (10) flips do not occur and the energy is conserved. 
Therefore 


f(é =&, for &1 —&) <0. (15) 
In particular, — 2/3 < — < 0 and — > 1 correspond to the states of oscillation and 
p ] 


circulation respectively. 
A few curious details of the motion will now be summarized. 
1. “Oscillatory” flips, for which the sense of @; is opposite to that of 6; , occur if 


0<£< 3/3 —3']'” = 0.56; “circulatory” flips occur if 0.56 < & < 1. The separating 
value corresponds to @; = 0. 

2. The maximum energy loss per flip corresponds to £ = 3°'”* = 0.57, with the 
parabola passing through the point of suspension. 

3. The “amplitude”’ of a flip can be defined as a = | sin [(6, — 6,)/2]). Its maximum 
a = 1, corresponds to & = 2°” = 0.71, 6, = 3/4, 6, = 2/4; its minimum, a = 0, 
occurs when é = 0, 6, = 6. = 7/2 and when é = 1, 6, = 6, = =. 


4. The entire history of the motion can now be described in terms of the function 
f(x), which is graphed and tabulated below, together with the corresponding values of 
6, and @, . 


f (x) 








To Te 





Fic. 1. Graph of the function f (2). 
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TABLE 1 


Function f(xz),O0 <2 <1 











M f(z) | 6; 62 

0.0 0.000 90°.0 90°.0 
0.1 0.079 95 .7 72 .8 
0.2 0.049 101 .5 55 .4 
0.3 —0.134 | 107 .5 37 .6 
0.4 —0.409 113 .6 19 .3 
0.5 —0.625 | 120 .0 0.0 
0.6 —0.608 | 126 .9 —20 .6 
0.7 | 0.271 134 .4 | —43 .3 
0.8 0.290 143 .1 —69 .4 
0.9 0.793 154 .2 —102 .5 
1.0 1.000 | 180 .0 —180 .0 





Some of the important properties of f(z) are listed below: 
1. f(x) is of class C’; i.e. it has continuous derivatives up to the second order in its 


entire domain — 2/3 < —E< o~, 
2. The stationary values of f(x) are given by min f(x) = — 0.655 = f(0.544), max 
f(x) = + 0.086 = f(0.130). 


3. The real zeros of f(a) are: 
zx, = 0, Xe = 0.236, x; = 0.751. 
1. These zeros, together with the boundary point x, = — 2/3 and the recursive 
relation 
f(Xn+0) = Less ; k=0,1,-:-, (16) 


define a monotone increasing sequence {z,};k = 0,1, --- , converging to x. = 1. This 
result follows from the facts that z, and 2, lie in the interval max f < x < 1, that in 
this interval the inverse function f-‘(x) is defined and is bounded by the inequality 
x < f ‘(x) < 1. The first nine terms of the sequence are: 


—0.667, 0.000, 0.236, 0.751, 0.791, 0.889, 0.899, 0.930, 0.935,---, (17) 


The value é lies in some interval J, , such that 


t; S28 << Biss if —2/3 <& <1, (18) 
kw e if 1<& <0. 


The three intervals k = 0,k = 1,k = © are the only stable states of ¢, in the sense that 
the system once in such a state will remain in that state forever. In particular, the 
two extreme states, k = 0, — 2/8 < § < O,andk = ~©,1< & < — are the classical 
oscillation and circulation respectively. The state k = 1,0 < & < 0.236 contains an 
infinite succession of oscillatory flips whose amplitude converges to zero; the motion 
converges to an oscillation with £ = 0, max @ = 2/2. In all other states a flip results in 
the jump Ak = — 2 in the step-function k(t). Consequently the particle executes 
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[k(0)/2] flips before descending to the state k = 0 if k(0) is even or to the state k = 1 
if k(0) is odd. The bracket above denotes the integral part of a number. The total number 
N of flips and the ultimate state k() of the system are hence given by the following 
table. 

TABLE 2 


History of the motion 


k(0) N k( «) Ultimate state 
even [k(0) /2] 0 oscillation 
odd co ] flipping 
2 0 circulation 
As an example, consider (0) = 0.777. From inspection of the sequence (17) it is 
seen that k(0) = 3. After one flip the system reaches the state k = 1, characterized 


by an infinite succession of flips, as indicated in Table 2. It is to be observed that if 
the support were rigid the pendulum would circulate, since ¢ exceeds the critical value 
£ = 2/3, (ec = 1), which separates the oscillatory and the circulatory states in the classical 


= °, = 


case. 
4. Summary. In the case of a simple pendulum supported by a flexible cord, when- 


ever the tension vanishes the particle passes from the circular to a parabolic trajectory. 
The energy loss occurring upon re-entry into the circle is given by the expression 
64 
3 
= 0, (x(1 — x) < 0), 


[x(1 — 2°)’, 0 <2 <1) 


where x and f(z) are the old and the new values of the energy, measured from the hori- 
zontal configuration @ = 7/2, and normalized by a divisor 3 mgl/2. The state of the 
pendulum is characterized by an integer k, such that 


x = a a : —— 
k= @, r> 1), 
where the sequence {z,};k = 0, 1, «-- , is monotone increasing, converging to 3 
and defined by 
f(Ti+s) = Leve k= 0,1 
Lo = —§, a, = 0, 1%. = 0.236, a, = 0.751. 


Here Zs, Fes and ZT; are the three real zeros of {(x), and x, is the lowest possible energy. 
The two extreme states k = 0 and k = © correspond respectively to pure oscillation 
and pure circulation. The state k = 1, which does not arise in the usual case of a rigid 
support, contains an infinite succession of parabola-circle transitions. The three states 
k = 0,k = 1 andk = © are the only stable states, in the sense that the system once 
in such a state will remain in that state forever. From an unstable state the system will 
ultimately descend to the state k = 0 if k(0) is even, or to the state k = 1 if k(0) is odd. 
In this process there occurs a finite number of parabola-circle transitions, this number 
being equal to [k(0)/2], the largest integer not exceeding k(0)/2. 
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NOTE ON A PAPER BY J.R.M. RADOK* 
By IAN N. SNEDDON (The University of Glasgow, Scotland) 


It is of interest to note that the complex variable solution of the equations of dynamic 
plane elasticity recently derived by Radokf can be arrived at by splitting the displace- 
ment field by the classical decomposition 


ay ,_ _av, a0 


pS — 9 . 
Ox Oy Ox oy 


and noting that the wave equations 
aVeo=¢ , aVip =" 
to which they lead have solutions of the form 
d(x — ct + iB,y), v(x — ct + 78.y) 


with 8? = 1 — c’/c} , 83 = 1 — c’/c} . The details are given in a paper published a few 
years ago**. Radok’s solution (4.4) is immediately derivable from my equation (8) by 
replacing my f’ and g’ by —¢/z and 7(1 + 83)~/(2u8,) respectively. 


AN ALTERNATE SOLUTION OF STEFAN’S PROBLEM{ 


BY 


ARTHUR L. RUOFF (Cornell University) 


The author, using a simple transformation, has solved two cases of Stefan’s problem. 
These problems have been solved earlier by Stefan [1] and by Neumann [1] by other 
methods. There are a number of problems governed by the parabolic differential equation 
with a moving boundary condition. Such a problem is the melting of a solid in which 
ease a finite heat sink exists at the position of the moving boundary. A similar problem 
exists in certain crystalline transformations. These problems, known as moving boundary 
problems or Stefan’s problems, also arise in studies of gravity drainage and seepage of 
oil from sand beds during pumping. 

The first problem considered here is the same as that studied by Stefan, being the 
ease of a solid initially at the melting point with one face brought instantaneously to 
some temperature above the melting temperature and held constant thereafter. The 
remaining sides are considered as insulated so that the heat flow is one dimensional. 

The second problem is the same as that discussed by Neumann, i.e., the moving 
boundary problem for a semi-infinite bar initially at a constant temperature below the 

*Received July 24, 1957. 

tJ. R. M. Radok, Quart. Appl. Math. XIV, 289 (1956). 

**I. N. Sneddon, Rend. Cire. Mat. Palermo, (77) 1, 57 (1952) 

tReceived January 31, 1957. 
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melting temperature with the leading face brought to some higher temperature instan- 
taneously. 

Solution of Case I. Consider a rod of ice at 32°F. Bring the face of the rod to some 
temperature u(0, ¢), the other faces being insulated. At a given time the equation 


Vu(x, t) du(a, t) 
Pula DL Sue, (1) 
Ox ot 
where 
a = k/C,p 


describes the temperature in the liquid region. The boundary conditions are, say, 


u(O, t) = 0, (2a) 
u(x(t), 0) = C, (2b) 
@) ae 
—k = = pAH ,x'(t), 
Ok nnets 


where C is the temperature of fusion, x(t) is the distance of the ice-water interface from 
the front face, p is the density, AH, is the heat of fusion and x'(t) is dx/dt or the rate 
at which the interface moves. 

This last equation may be written as 


Ou a 
—-— = — 2 (2), (2c) 
OF Newsti a ry F 


where A = AH,/C, . 
Transforming Eqs. (1) and (2) in terms of 


y(2, t) = nae (3) 
we obtain the equation 
i + 2y = = 0 (4) 
with the boundary conditions 
u(0) = 0, (5a) 
u(y()) = C, (5b) 
= = 2A[y(t) + 2ty |. (5e 
dy Sash 
By direct quadrature Eq. (4) yields 
o Be~*" (6) 


and 
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Applying (5a) and (5b) to the latter gives 


y(t) 
C= 8 [ e” dy. (8) 
Jo 
Hence y(t) = y, a constant, since C is a constant. This immediately gives the time de- 
pendence of the position of the moving boundary, 
x(t) = 2y(at)"”?. (9) 
Equation (5c) now becomes 
—Be~*” = 2Ay (10a) 
and Eq. (8) becomes 
7 . 
C= e | e” dy. (10b) 
0 


The equations (10a) and (10b) can be readily solved for y and @ by first eliminating 3 and 
solving the resulting equation for y since the distribution function and the error function 
are both tabulated. 

Finally 


au 


u=8 | e" dy (11) 


0 


subject to 0 < y < y, a condition which requires that the solution hold in the liquid 
region only. 

Solution of Case II. Dividing the bar into region I which is the fluid region and 
region II the solid region, we obtain the following conditions on the rate of heat flow: 





du; ad Ouy 

“Tox* at’ (12) 
u(0, f) = 0, (12a) 
u,[x(t), t] = C, (12b) 
Q os “i - 
i, <p aelh ~ k te, (12¢) 

Ox Ox 
un _ Jun (1 
Oy] eC? (13) 
Ur; [x(0), t] = ie (13a) 
t3(0, ft) = az, 0) = T,, (13b) 
x(0) = 0. (13¢) 


For convenience it is assumed that the density does not change. A further assumption 
is that equilibrium conditions are satisfied at the interface as regards temperature and 
state of aggregation; since the process of melting is a rate process it is expected that the 
above approximation would not be valid for extremely rapid heating rates. The problem 
of recrystallization of metals is often complicated in such a way. For example, the trans- 
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formation from the face-centered cubic form in steel to the body-centered form is a 
rather slow process so that the condition (12c) is further complicated by the addition 
of some time-dependent factor on the expression pAH,x'(t). There are, however, crystal- 
line transformations which fall into the category described by Eqs. (12 and 13), these 
being the metals in which pure martensitic-type transformations occur. 

In order to solve the given equation let 


z=2/t (14) 
Equation (12) now becomes 
du 2duy is 
a ge? * Ode ~ ie 
while (12a, b and c) become 
u(0) = 0, (15a) 
ul2)] = C, (15b) 
| 
—k, is = pAH (tz + 2/2) — ky os , (15e) 
dz GZ \enuce 
Similarly Eq. (13) yields 
ox a re su - 0 (16) 
with the conditions 
urrle()] = C, (16a) 
“deo = T,. (16b) 


We leave the condition 2(0) unspecified for the time being. Treating du,/dz as the 
dependent variable Eq. (15) can be directly integrated, yielding 


du, 3 —2* (17 
= B exp ——-: 7) 
dz P day, 
tepeated integration gives 
C= | du; = B; | exp —— de. (18) 
Jo : fay 
This implies that z(t) = y, a constant. This then gives for the temperature in the 
i 
liquid region 
Uy = By | exp —— dz, (19) 
P 4a; 


where 8; is yet to be evaluated. Moreover, it gives the position of the moving boundary 
as a function of time 


alt . yt'”?, (20) 


where y is yet to be evaluated. 
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Integration of Eq. (16) is carried out in the same way, yielding 


2 








duyy —2z 
a exp 21 
dz Bir J dary ( ) 
and 
To ao —_—s 
/ dus; = Bis | exp —— dz (22) 
Cc J2(t) 4a, 
where z(t) = y. 
Equation (22) may be rewritten in the form 
0 er 4 — 
T —C = By [ exp —— dz — By [ oa —— a. (23) 
Jo fay, Jo dary 


The heat balance condition (15c) at the moving boundary has not yet been used. 


2 


—k,B; exp — = eh: _ kB exp ean (24) 
4a; 2 dary 


Equations 18, 23 and 24 are sufficient for obtaining 8,, 8, and y. From Eq. 18 we obtain 


a — -1 
6, = dl | exp rong ae| (25) 


and Eq. (23) yields 


no a? 9 -1 
By = (T. — of | exp —— dz — [ exp -—— ie | . (26) 
7 dary J ( tay; 


0 
Substitutions of (25) and (26) into (24) gives an equation involving y only. This can be 
solved by trial and error for y. Then 8, and 8;; are readily available from Eqs. 25 and 26. 

Finally, the position of the moving boundary is given by Eq. 20 and the expression 
for the temperature is given for regions I and II, respectively, by 


— 6: | exp — dz, (27) 
0 da; 
where z < 2(t) = ¥, 
and, 
Ui; = To Br | exp — dz (28) 
Jez dary 
with z > z(t) = y. 
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BOOK REVIEWS 


(Continued from p. 154) 


Cours de Mathématiques. By J. Bass, with a preface by G. Darmois. Masson, Paris, 

1956. xi + 916 pp. $24.55. 

The comprehensive text is based on the author’s courses at the Ecole Nationale Supérieure de 
]’Aéronautique and the Ecole Nationale Supérieure des Mines in Paris. These courses are designed to 
acquaint students of engineering with the mathematical tools of the modern engineer. 

As a detailed survey of the contents of the volume would far exceed the space available for this 
review, the following condensation of the table of contents will have to serve as an indication of scope: 
Linear Algebra (vector spaces, matrices, Hermitian space, tensors); Integrals (sets and integrals, defi- 
nite integrals, reduction of integrals, numerical quadrature, improper integrals); Functions Defined by 
Series or Integrals (Fourier series, orthogonal functions, Fourier integral); Curves, Line Integrals 
(differential geometry of curves, line integrals, application of line integrals to mechanical quadrature); 
Surfaces, Multiple Integrals (double integrals, change of variables, elements of differential geometry 
of surfaces, multiple integrals, surface integrals, integral formulas of vector analysis, improper multiple 
integrals, gamma and beta functions, elements of operational calculus); Analytic Functions (derivative 
of a function of a complex variable, elementary analytic functions, integration, series, theorem of resi- 
dues); Ordinary Differential Equations (general properties, linear differential systems, linear differential 
equations of the second order, Bessel functions, numerical integration); Partial Differential Equations 
(linear first-order equations, non-linear first-order equations, linear second-order equations, examples 
(vibrating strings and membranes, sound waves, heat conduction), Laplace equation and Newtonian 
potential); Appendix (calculus of variations, nomography). 

The exposition is clear and easy to follow; results are always stated with useful generality and rigor 


even if, for the sake of brevity, they have not been derived in full generality. 
W. PRAGER 


Science and information theory. By Leon Brillouin. Academic Press, New York, 1956. 

xvii + 320 pp. $6.80. 

This book is a series of essays on topics in the theory of information, in the physics of thermody- 
namics and fluctuations, and on relations, established or surmised, between these domains of discourse. 
The essays themselves are based on lectures given by the author to engineers of the International 
Business Machines Corporation, and upon other lectures and papers by the author. 

About one-third of the book, Chapters 1 through 8, is devoted to aspects of the theory of informa- 
tion. The coverage is restricted to the introductory material needed for the author’s later discussion, 
and to some special topics. A comprehensive introduction to the theory is not attempted. 

The remainder of the book is devoted to the author’s main thesis—that the entropy of physics and 
Shannon’s measure of information are one and the same. His treatment of this thesis is that of the physi- 
cist; it is in the same spirit as several briefer studies of the subject by others, studies to which the author 
makes generally adequate references, 

In more detail: The book’s first four chapters (50 pages) present in the author’s words some of 
the elements of Shannon’s first paper on his theory of communication. (“A Mathematical Theory of 
Communication’ Bell System Technical Journal, v.27, pp. 379-423, pp. 623-656.) 

Chapters 5, 6, and 7 (28 pp.) treat the minimal redundancy codes of Shannon and Fano, and the 
error correcting codes of R. W. Hamming. The author exhibits some interesting alternatives of his own 
to these latter. Shannon’s formula for the capacity of the binary channel with symmetric, independ- 
ent, and memoryless noise is introduced, but the general concept of channel capacity is not defined. 

Chapter 8 (36 pp.) discusses some topics in Fourier analysis which are useful for studying the sta- 
tistical aspects of electrical communication: the uncertainty relation between time and frequency, the 
so-called “sampling theorem” for a function whose Fourier transform vanishes except over an interval, 
(electrical) filters, the autocorrelation of a function (introduced as a time average, not an ensemble 
average), and the Wiener-Khintchine representation thereof. 


1958] BOOK REVIEWS 203 


on 
Qo 


Chapters 9, 10 and 11 (38 pp.) present a summary of phenomenological thermodynamics, and dis- 
cussions, both phenomenological and in the idiom of Gibbs, of Brownian motion and of thermal fluc- 
tuations in simple mechanical and electrical systems. 

With Chapter 12 (11 pp.), ‘The Negentropy Principle of Information,” begins the author’s 
discussion of his main theme: the principle to which he has given the name just quoted. Chapter 12 
enunciates the principle, with illustrations, as a generalization of the Second Law of thermodynamics. 
Chapter 13 (24 pp.) applies the principle to the special problem raised by Maxwell’s Demon, who, by 
defying the Second Law for almost ninety years, has outraged three generations of physicists. The 
author exercises this demon by invoking his negentropy principle, in this way going somewhat beyond 
the more common observation that the demon could not function, at least anthropomorphously, if 
the local radiation field were in thermal equilibrium with the ensemble he was upsetting. The argu- 
ments in these two chapters are a mixture of phenomenological thermodynamics and quantum kinetic 
theory. The reviewer interprets these chapters as the author’s proof that his principle is valid, a proof 
obtained by showing that the principle is necessary in order to preserve the Second Law, in particular 
against daemonic attacks. 

Chapters 14, 15, and 16 (62 pp.) apply the negentropy principle to problems of measurement in 
many physical contexts. Again, the methods are a mixture of phenomenology and appeals to the princi- 
ple of equipartition in kinetic theory. 

Chapter 17 (16 pp.) applies the negentropy principle specifically to the physical problem of com- 
municating electrically in the presence of thermal fluctuations. The argument appeals to equipartion; 
it results in the formula of Shannon for the capacity of a channel with additive Gaussian noise. A defi- 
nition of channel capacity adequate for this particular case is enunciated. 

Chapter 18 (8 pp.) treats problems, troublesome in the light of the author’s principle when it is 
interpreted as a conservation principle, which arise because printed information can be duplicated and 
somehow multiplied by dissemination. These problems go far beyond the scope of any mathematically 
formulated theory of information, but seem to yield to the author’s phenomenological methods. 

Chapter 19 (18 pp.) treats some topics in computing, and Chapter 20 (14 pp.) several unsolved 
problems, some of which also go far beyond any presently formulated mathematical theory. 

Just as there are, broadly, three views of thermodynamics, the phenomenological, that of kinetic 
theory, and that of Gibbs—so there can be three views of the author’s negentropy principle of informa- 
tion. It is difficult for a mathematician to catch the author in a fully satisfactory statement of any of 
these views. The reviewer will attempt some clarifications. 

First, a word about signs. Shannon regards — = p; log pi, quite naturally, as a measure of “un- 
certainty.” When this “uncertainty” is eliminated or reduced, “information” results. Both quantities 
are treated as being intrinsically positive, just as a solvent person who is not an accountant naturally 
regards both his bank balance, and the pocket money he has just withdrawn from that balance, as 
intrinsically positive. Brillouin follows this usage with good consistency, and does not display the con- 
fusion of some of his successors on the matter. To avoid here any such confusion, interpret “information” 
as “information obtained as the result of a physical change.” When this quantity is positive it represents 
a withdrawal from a balance of “uncertainty.” Therefore “negative information” is, from the point 
of view of an accountant, “uncertainty.” 

Phenomenologically, the negentropy principle of information is a postulate, a conservation postu- 
late, in which ‘“‘negative information” (i.e. “uncertainty” in the reviewer's words above) is introduced 
as a balance to entropy. In a reversible transformation of a closed system the sum of “negative infor- 
mation” plus entropy then remains constant. At this point, “negative information” appears, like po- 
tential energy and the neutrino, as a concept physically necessary to preserve a conservation law. In 
this way the conservation postulate can be regarded as defining “physical information” (reviewer's 
term introduced for convenience). 

Brillouin chooses to extend this formulation, as is usually done with the restricted form of the 
Second Law, to an inequality applicable to irreversible processes: he asserts that, in a closed system, 
the change in entropy-minus-information (entropy plus negative information) is always non-negative. 
This is his enunciation of the negentropy principle as a generalization of Carnot’s principle (the Second 
Law). 

Since physical information is defined above by the conservation law only as a quantity whose 
negative balances entropy in a reversible change, the author must find some other, independently de- 
fined, form of “information” to insert in his generalization of the Second Law. For this, phenomeno- 
logical concepts do not suffice, and he turns essentially to kinetic theory. He chooses for his information 
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the quantity « In(Po/P:) where (page 152) Po is the “. . . number of possible cases” existing a priori, 
and P; the number of cases which are possible after the physical change has taken place. Both a priori 
and a posteriori, the cases possible are postulated to be all equally likely. The numerical factor « is, 
quite naturally, Boltzmann’s constant (dimensions: energy per temperature). 

After this essentially kinetic statement of the negentropy principle, three arguments, or kinds of 
argument, intertwine. The first is the phenomenological argument of Szilard that the conservation 
postulated above, (i.e., that physical information is a balance to entropy) is necessary to preserve the 
Second Law. The second is an argument from kinetic theory defending «In(Po/P;) as that which bal- 
ances entropy in this conservation: by this then the author identifies what above was called physical 
information with concepts from kinetic theory. Finally, there is much use of the language of the theory 
of communication to associate « In(Po/P1) with concepts found in that theory. 

A complete formulation of the author’s negentropy principle in terms appropriate to kinetic theory 
would be mathematically quite interesting. The subject is too complex for discussion in a book review, 
but it is evident that the author has missed some of the interesting features. He has, for example, failed 
to observe that in kinetic theory entropy appears most naturally not as a function of state variables 
nor as a functional of the ensemble, but as a random variable—in fact, that random variable whose 
expectation is Shannon’s average conditional entropy of the complexion given the macroscopic state. 
Once one recognizes this, Shannon’s theorem on the asymptotic approximate constancy of that variable 
applies to large systems and can be used to justify the equiprobability that has been assumed by Bril- 
louin, and by all physicists since Boltzmann. 

The author omits any discussion of the negentropy principle in terms appropriate to Gibbs’ view 
of thermodynamics. In such terms, the negentropy principle is simply this: in a canonical or grand 
canonical ensemble of (e.g., quantum) physical systems Shannon’s functional — Zp; log p;, in which p; 
is the probability that the 7 state be occupied, is, apart at most from a linear and monotone change of 
scale, equal to the thermodynamic entropy of that ensemble. In this form, the “principle” is tauto- 
logous, and it is false when applied to any but the canonical or grand canonical ensembles. From the 
Gibbs point of view, then, the interesting physical fact is not the negentropy principle but the universal 
usefulness of these canonical ensembles. 

This book suffers from hasty editing; many of its chapters are clearly derived directly from lecture 
notes, they overlap in places and often use with brief, or no, explanation ideas introduced more care- 
fully later. More seriously, the author has not defined his intended audience: a wide range of topics 
is covered, each in some detail and often with some kind of sophistication, but almost no topic is treated 
adequately for a beginner therein. 

The author is generous with his references. He does, however, overlook a particularly precise and 
satisfying treatment of Szilard’s argument that negative physical information (as the term is used in 
this review) balances entropy. This is the paper ‘“The Well Informed Heat Engine’ by Richard C. 
Raymond (American Journal of Physics, v.19, pp.109-112, 1951; compare its title with Section 6 of 
Chapter 13). 

As it appears in this book the substance “information” has several forms, some of them quite 
anthropocentric. In the typical argument, more than one form is displayed. There results an impression 
of elusiveness and the reader is often challenged to be sure that the author has indeed limited all appli- 
cations of his principle to strictly closed systems. In contrast, the less ambitious treatment by Raymond, 
just mentioned, is explicit and precise on this point. 

To a mathematician, it is clear that the author is a physicist. He exhibits no interest in logical com- 
pleteness—for example, Shannon’s full definition of channel capacity, and his fundamental theorem 
concerning it, are never mentioned. The author depends throughout upon phenomenological arguments 
in terms of ‘‘gedanken-experimente,”’ and upon the most naive formulations of the mathematical con- 
cept of information. Little effort is made to distinguish between random variables and quantities which 
define their distributions. Sample spaces, when they are mentioned at all, are specified only in the loosest 
terms; this latter despite the author’s own admission that “information” depends explicitly upon the 
sample space i.e., upon the ensemble of a priori possibilities. 

This reviewer feels that the author has a valid message and that a negentropy principle of infor- 
mation could be made mathematically precise. He regrets that, in failing to make the principle precise, 
Brillouin’s book has added the considerable weight of its author’s support to an already flourishing 
mystical school of information theory. 


B. McMILLAn 
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Quantum chemistry—an introduction. By Walter Kauzmann. Academic Press, New York, 
1957. xii + 744 pp. $12.00. 


The outstanding feature of this new introductory text on chemical quantum mechanics is its ex- 
cellent mathematical section, which is much more extensive than in most books of this sort. After pre- 
senting a certain amount of necessary background material, the author goes into a detailed treatment 
of vibrating systems, starting with the vibrating string, and continuing on to cover membranes of 
various shapes, the uniformly flooded planet, and three dimensional continua. In the course of this 
discussion all the basic mathematics of quantum mechanics, including perturbation theory and the 
variational method, is developed in terms of physical situations relatively familiar to the beginning 
student. 

The basic principles of quantum mechanics are presented next in a purely postulational manner 
which, while of course unassailable from a logical standpoint, is unlikely to appeal to most readers 
encountering this material for the first time. Following this there are two large sections on atomic and 
molecular quantum mechanics, well done in a more or less conventional fashion. A regretable omission 
here is the absence of any presentation of group theory, although symmetry considerations are ex- 
tensively used. 

The final part of the book is devoted to non-stationary situations, in particular the interaction of 
atoms and molecules with electromagnetic radiation, discussed in considerable detail both from clas- 
sical and quantum viewpoints. Here there is a good but somewhat lengthy treatment of optical rotation, 
given apparently at the expense of equally important topics such as the Raman effect, which is not 
mentioned at all. 

The exercises scattered throughout the text deserve especial mention. They are entirely original, 
workable without being trivial, and well chosen to clarify the accompanying material. 

On the whole, in spite of the few criticisms offered in this review, this is probably the best beginning 


text on quantum chemistry currently available. 


I 


STEPHEN PRAGER 


Viscous flow theory IJ—turbulent flow. By Shih-I Pai. D. Van Nostrand Co., Inc., New 


Jersey, 1957. xi + 277 pp. $6.75. 

Turbulence represents a field of fluid mechanics which is both very important and very difficult. 
rhe literature on the theory of turbulent flows ranges therefore from subtle mathematical studies to 
the crudest approximations. Experimental work on turbulence covers a similar range. To penetrate from 
the outside into the field of turbulence research through reading of original papers is difficult and the 
non-specialist or student trying to do this must find the literature very confusing indeed. 

\ text which covers the elements of turbulence research and the applicable results so far obtained 
is thus very much needed. According to the table of contents of Dr. Pai’s book it seems to fill at least 
part of this need. After closer study however I came to the conclusion that this is not the case and that 
the book in its present form will add to the confusion rather than help to eliminate it. 

The fundamental shortcoming in the book is the lack of a personal point of view of the author and 
of a definite scientific standard. Work with high mathematical standards and little physical insight, 
vork with much physical insight and less rigor, and crude ad hoc computations to check in one way or 
other experimental data are compiled indiscriminately. Thus in the chapters on the statistical theory 
of turbulence the author insists on using Lebesgue instead of Riemann integrals, and worries about 

ergodic theory. In computing shear flow problems on the other hand, computations which amount 
to little more than interpolation formulae for a specific set of experiments are presented as ‘‘theories”’ 


(e.g. Laufer’s measurements on pipe and channel flow). 

Reading the book I recognized chapter by chapter the impact of different personalities and schools 
of turbulence research which the author has neither assimilated into his own way of thinking or even 
style of writing, nor put into the proper context. The use of experimental results is equally unsatis- 
factory because experiments are used only to “verify” computations, and sometimes single experi- 
mental curves are taken out of their original context, thus defeating the purpose of the original work. 
Measurements on the turbulent boundary layer at high Mach numbers are cited in one sentence: 
“Most of the experimental values for the skin friction coefficient agree within about 10% with the the- 
oretical values of van Driest, Wilson and Rubesin et al.’’ Not long ago there existed over twenty dif- 
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ferent theories for the turbulent boundary layer at high Mach number and the choice of the one ‘‘most 
likely to succeed”’ is entirely based on a few sets of first rate experiments like Coles’ and Chapman and 
Kester’s. The author does not seem to realize that good theory like a good experiment can stand on its 
own feet and that agreement between the two does not improve a bad theory or a poor set of experiments. 

A number of misleading and meaningless statements which appear throughout the text can be 
traced to the same basic cause, i.e. the failure to assimilate opinions of others within a unifying point 
of view; e.g. (p. 20): ‘“‘The law of heat loss (of a hot wire at high speeds) is non-linear while that for low 
speed flow, King’s law, is linear.” Linear in what quantity?—(p. 109): ‘“The fact is that neither the tur- 
bulent Prandtl number nor the laminar Prandtl number are unity but rather about 0.75 for air.” Js 
the turbulent Prandtl number 0.75?—(p. 165): “The joint-probability distribution of fluctuating 
velocities at two points and the probability distributions of the derivatives of the velocity components 
are not normal, due to the fact that these quantities satisfy the Navier-Stokes equation; hence they 
are not entirely random in nature.’’ Cannot statistics be applied to any mechanical syst m? This list 
could be extended. Every one of these statements can be interpreted correctly by somebody familiar 
with the subject matter. For an outsider loose sentences like these must be puzzling and for a beginner 
utterly confusing. Here a good deal of editing would have helped. 

As far as the selection of subject matter is concerned the overall choice of the chapters is in my 
opinion quite good. However everywhere one finds detailed manipulation of equations rather than 
critical and constructive explanations of why certain assumptions make sense and others not. Thus 
four different approaches to turbulent boundary layer theory at high Mach number are given in great 
and cumbersome detail; on the other hand one of the most fruitful concepts of recent years, Clauser’s 
“equilibrium layer,’ is barely mentioned. I also miss Donaldson’s simple and intuitive approach to 
Mach number effects on turbulent skin friction. In general the simple and general laws which would 
be most helpful to the reader are not explained: e.g. the law of the wall, the linear stress-distance re- 
lation in pipe and channel, the similarity solution for free turbulence. All of these are there but they 
are smothered under lengthy manipulations of equations. Here too, I believe, one notices again that the 
book is assembled rather than written. 

Summarizing, I regret that I cannot recommend the book as it stands since it will tend to give the 
reader a distorted and often misleading impression of what turbulence research is and what has been 
achieved. To read up on turbulence the student still has to be referred to Goldstein’s ‘‘Modern De- 
velopments” and to some survey papers and original work until he is ready to read Batchelor’s and 


Tow isend’s recen m 1ograpns. 
H. W. Lik PMANN 


Applied probability. Proceedings of Symposia in Applied Mathematics. Vol. VII. Mc- 
Graw-Hill Book Co., Inc., New York, Toronto, London, 1957. v + 104 pp. $5.00. 


This volume contains nine papers presented at the seventh symposium in Applied Mathematics 
of the American Mathematical Society. The intention of the symposium was to explore three principal 
subjects, the theory of diffusion, the theory of turbulence and probability in classical and modern 
physics. No attempt is made to review these very extensive subjects although a few papers review 
special topics and others deal with problems of more recent origin than the above titles might suggest. 
The papers describe work of the high quality that would be expected from the following very impres- 
sive list of authors: Paul Lévy, J. L. Doob, William Feller, Eberhard Hopf, Guido Miinch, G. K. 


Batchelor, Mark Kac, 8. M. Ulam, and B. O. Koopman. 
G. F. NEwELL 


Relaxation methods in theoretical physics, Vol. II. By R. V. Southwell. Oxford University 
Press, New York, 1956. 274 pp. $8.80. 


This book presents the latest developments in the solution of boundary value problems in partial 
differential equations by the numerical, hand computation scheme originated in the 1930’s by the 
author. The book is a second volume on this subject, the first volume of which appeared in 1946. These 
two volumes, together with the book ‘‘Relaxation Methods In Engineering Science” which appeared 
in 1940, present essentially all that has been accomplished in development of the method up to 1955. 

The method of solution basically is to convert the differential equation to a finite difference net 
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and to assume a numerical value as the solution at each net point. These incorrect values are then 
changed at one point at a time in such a way that the successive values approach the correct ones. 

In the present volume, as in the past ones, the methods of carrying out the solution in detail are 
discussed by specific examples. In large part the illustrations are taken from elasticity; stresses and 
deflections of plates and solids of revolution, vibrations of plates and membranes, and elastic stability. 
A discussion of electromagnetic oscillations forms an obvious application of some of the same methods. 
As an example the lower modes of oscillation of a Klystron tube are computed. 

Chapters XIII and XIV (Vol. II contains chapters VII to XV) deal with particularly difficult 
non-linear problems; large deflections of plates, the flow of viscous fluids and plastic strain. 

The final chapter makes brief mention of a number of recent, only partially developed ideas. Some 
three dimensional problems are solved and methods of applying relaxation techniques to parabolic and 
hyperbolic type equations by converting them to elliptic equations of higher order are noted. 

The illustrative problems and many of the ingenious tricks of solution are taken with acknowledge- 
ment from a number of workers. Many of these are former students of the author and have collaborated 
with him in making the present as well as the former volumes possible. 

For many problem solvers who have become accustomed to think in terms of the large scale com- 
putation machine, the relaxation method may appear like a good idea that arrived a century too late. 
From the nature of the problems solved, it is clear that considerable time must be spent with a desk 
computer to find them. However, there are many people who for years will not have a readily available 
large scale computing machine. Furthermore, some of the problems solved have a complexity calling 
for judgment during solution which is not yet available in a computing machine. At present, of course, 
computing machines sometimes use an identical finite difference net with an iterative method of solu- 
tion. If a future machine is ever built which can scan a field of numbers as much faster than it can 
compute as is true of the human computor, the relaxation method will emerge as a most significant 
development and the present volume will be one of the classics. 


H. W. Emmons 


Digital computer programming. By D. D. McCracken. John Wiley & Sons, Inc., New 
York, 1957. vi + 253 pp. $7.75. 


A first book on any subject is always welcome, even if it did not fulfill its purpose as admirably 
as does the book under review. 

It is elementary, to be sure. A book is still to be written (or, at least, published) on modern high 
speed computing methods applied to advanced problems in numerical analysis. But there has not even 
been available before a systematic introduction to the art of coding for an automatic stored program 
computer. The computer on which we learn to solve problems is TYDAC, a mythical computer whose 
approximate counterpart in real life the experts will easily be able to guess: it is a decimal machine 
with 2000 words high speed storage, index registers and four tapes. The discussion includes number 
systems, fixed and floating decimal point methods, address computation with and without index regis- 
ters, loops, subroutines, flow charts, input-output methods, and all too short an anatomical discursion 
into interpreters and compilers. The explanations are uniformly lucid, the presentation is systematic 
and purposeful, the layout excellent. 

As a textbook for a first programming course and as a compendium to numerical analysis seminars 
the book is heartily recommended. 

W. F. FREmERGER 


Automatic digital computers. By M. V. Wilkes. John Wiley & Sons, Inc., New York, 
1956. x + 305 pp. $7.00. 


The contents of this volume provide a general introduction to the principles underlying the design 
and use of digital computers. After a historical introduction the author, who writes with authority of 
many years pioneering effort in the field, presents the principles of logical design and of program con- 
struction. The latter are illustrated by means of the order-code of the Cambridge EDSAC. Storage, 
electronic switching and computing circuits, design and operation of digital computer, both relay and 
electronic, are discussed in considerable detail. 
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The book is by no means easy to read, but the expert will find in it a wealth of information other- 
wise hard to come by. The solution of particular problems is not discussed—the application of com- 
puters to the solution of advanced problems in numerical analysis is, by self-imposed limitation, outside 
the scope of the work which is thus of relatively little interest to mathematicians. Its value lies in giving 
an up-to-date picture, at the time of writing, of the state of computer design on both sides of the 


Atlantic. 
W. F. FREIBERGER 




















